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Chapter 1 

The ecology of black holes in star 
clusters 

S. F. Portegies Zwart 
University of Amsterdam 



1.1 Introduction 

In this lecture we investigate the formation and evolution of black holes in 
star clusters. 

The star clusters we investigate are generally rich, containing more 
than 10 4 stars, and with a density exceeding 10 4 stars/pc 3 . Among these 
are young populous clusters, globular cluster and the nuclei of galaxies. 

Under usual circumstances black holes are formed from stars with a 
zero-age main-sequence (ZAMS) mass of at least 25-30 M Q (Maeder 1992; 
Portegies Zwart, Verbunt & Ergma 1997; Ergma & van den Heuvel 1998). 
These stars live less than about 10 Myr, after which a supernova results in 
a black hole in the range of 5-20 M (Fryer & Kalogera, 2001). The rest of 
the mass is lost from the star in the windy phase proceeding the supernova 
or in the explosion itself (Heger 2003). For a Scalo (1986) or Kroupa & 
Weidner (2003) initial mass functions (IMF), variants of Salpeter (1955), 
one in 2300-3500 stars collapse to a black hole. So black holes are rare, but 
still, the Milky-way Galaxy is populated by 30 to 40 million black holes. In 
the remainder we will refer to these relatively common type of black holes 
as stellar mass black holes or simply as BH. 

Supermassive black holes (SMBH) are, with about one per galaxy, 
among the rarest single objects in the universe. They have masses of about 
10 6 to 10 10 M© which also makes them among the most massive single ob- 
jects in the Universe. The best evidence for the existence of a supermassive 
black hole comes from the orbit of the star S2 which is in a 15.2 year orbit 
around an unseen object (Schodel et al. 2002). The mass of this black 
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hole is 3.7± 1.5 x 10 6 M Q , which puts it directly at the bottom of the mass 
scale for SMBH's (see also fig. lTT7) . 

The gap in mass between stellar mass black holes and supermassive 
black holes may be bridged by a third type, often referred to as intermediate 
mass black holes (IMBH). 

In this chapter our main interest is the formation and evolution of 
these black hole families in star clusters. Also the possible evolutionary 
link between stellar mass black holes, via intermediate mass black holes to 
supermassive black holes will be addressed. We mainly focus, however, on 
the ecology of star clusters. The term star cluster ecology was introduced 
in 1992 by Douglas Heggie to illustrate the complicated interplay between 
stellar evolution and stellar dynamics, which by their mutual interactions 
has similarities with biological systems. 

In the here studied ecology we mainly address the gravitational dy- 
namics, stellar evolution, binary evolution, external influence and how these 
seemingly separate effects work together on the star cluster, much in the 
same way as in an organism. 

1.1.1 Setting the stage 

The main objects of our study are clusters of stars. There are a number 
of distinct families of star clusters, with one common characteristic: a 
star cluster is a self gravitating system of stars, all of which have about 
the same age. We make the distinction between various types, among 
which are: open cluster like Pleiads (see fig. ll.il) . young dense cluster like 
R 136 in the 30 Doradus region of the large Magellanic cloud or the Arches 
cluster (see fig. ll.2[) . globular cluster like M15 (see fig. ll.3fl . If galactic nuclei 
contain a population of stars with similar ages we could include these in the 
definition. The nucleus of the Andromeda galaxy, depicted in Figur eTHI 
is an example, though the stellar population in its nucleus has probably a 
large spread in ages. 

Figur eTTTl shows the ~ 115Myr old star cluster Pleiads 1 , at a distance 
of about 135 pc (Pinfield et al. 1998; Raboud & Mermilliod 1998; Bou- 
vier et al 1998). This cluster contains about 2000 stars, half of which are 
contained in a sphere with a radius of about 8pc. 

Figur eTl"2l shows Hubble space telescope (HST) images of two well 
known young dense star clusters, Arches (to the left) and NGC 2070 (right) 
in the Large Magellanic Cloud. These clusters are the prototypical exam- 
ples of young dense star cluster (YoDeCs), which are young ( ^S lOMyr), 

1 According to the Greek mythology, one day the great hunter Orion saw the Pleiads as 
they walked through the Boeotian countryside, and fancied them. He pursued them for 
seven years, until Zeus answered their prayers for delivery and transformed them into 
birds (doves or pigeons) placing them among the stars. Later, when Orion was killed 
(many conflicting stories as to how), he was placed in the heavens behind the Pleiads, 
immortalizing the chase. 



Introduction 




Figure 1.1. The Pleiads open star cluster (image by David Malin) contains a 
few thousand stars in a volume of a several parsec cubed. 





Figure 1.2. Left: Arches star cluster (HST image by Figer et al. 1999). At a 
projected distance of ~ 30 pc from the Galactic center this cluster is strongly in- 
fluenced by external forces. The cluster on the image still looks quite symmetric 
since we only observe the inner part. 

Right: The young dense star cluster R136 (NGC 2070, HST image by N. Wal- 
born)in the 30 Doradus region of the large Magellanic cloud. 
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massive ~ 10 4 M and dense n £ 10 5 stars/pc 3 . The two clusters show also 
that the class of young dense star clusters hosts two families of star clus- 
ters; the strongly tidally perturbed cluster (like Arches) and the isolated 
cluster (like R136). 




Figure 1.3. The globular cluster M15 (NGC 7078, HST image by 
Guhathakurta) . 

Figur eTTol give a Hubble's view of the core collapsed globular cluster 
M 15 (Guhathakurta 1996). It is one of the old globular star clusters in the 
Milky-way's halo, contains about a million stars, and the central density 
is of the order of n ^ 10 5 stars/pc 3 . No core has been measured in the 
density profile, and therefore this cluster is identified as a collapsed globular 
(Harris 1996, http: //phy sun. physics .mcmaster . ca/Globular .html). 

In figure rni we show an image of the Andromeda Galaxy 2 A total of 
693 candidate globular clusters were found in recent 2MASS MIR obser- 
vations (Galleti 2004). It is still unclear why M31 has more than trice as 
many globular clusters than the Milky-way Galaxy. Note also that there 
are no YoDeCs observed in M31, which is also somewhat puzzling, as the 
Milky-way Galaxy has at least four. 

In figure lTTBI we place a few well known star clusters in retrospect. 
All cluster images are on the same scale. It is interesting to see is how 
large the globular cluster M80 is compared to some depicted open clusters 

2 Andromeda was the Ethiopian princess, whom Perseus rescued and married. She 
became queen of Mycenae and, after her death, a constellation. The galaxy was later 
named after the constellation. 
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Figure 1.4. The Andromeda galaxy (M31, image by Jason Ware). 

(Pleiads) and compared to the known Galactic YoDeCs (Arches, Quintuple, 
NGC3603 and Westerlund 1). 

There appears to be a clear relation between the number of globular 
clusters, like M15 (see fjg. 11.30 . and the Hubble type of the host galaxy. 
This relation is expressed in the specific number of globular clusters Sn- 
For open clusters and YoDeCs no such relation is known, but young star 
clusters are particularly abundant in starburst and interacting galaxies like 
the Antennae and M82. 

Table lT"2l lists the space densities and specific numbers of globular clus- 
ters Sn per M v = — 15 magnitude (van den Bergh 1984) for various Hubble 
types of galaxies. The values given for Sn in Table ll.2l are corrected for in- 
ternal absorption; the absorbed component is estimated from observations 
in the far infrared. 

s N = iV GC 10 - 4(M " +15) . (1.1) 

Here NqC is the total number of globular clusters in the galaxy under 
consideration. The estimated number density of globular clusters in the 
local Universe is 

(p G c = 8.4 h 3 Mpc- 3 (1.2) 
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Trapezium 



•t4jfc : _ Westerlund L 
Arches 
R136 
Quintuplet 



Figure 1.5. Star clusters in the correct physical scale. The Galaxy fig. ll.4l was 
left out for practical reasons. Images are from: M80 (F. Ferraro); Trapezium 
cluster (HST image by . Bally, D. Devine, R. Sutherland, and D. Johnson); 
Westerlund 1 (2MASS image by S. van Dyk); Arches (HST image by Figer, 
1999); R136 (HST image by N. Walborn); Quintuplet (HST image by Figer 
1999); MGG-11 (VLT image by McCrady et al 2003); Pleiads (AAT image by D. 
Malin). 



(where h = Ho/100 kms Mpc 1 ), slightly smaller than the result re- 
ported by Phinney (1991). 

Figur eTl"r)l illustrates the concentration of clusters, expressed in the 
structural parameter Wo (King 1966), which ranges from about 1 for very 
shallow clusters to about 12 for very concentrated clusters. The hgure 
shows the images of three clusters with quite different concentration ranging 
from Wq — 6 (Omega Centauri) to Wq ^ 12 for the core collapsed globular 
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Table 1.1. Overview of selected parameters for young dense clusters (Massey & 
Hunter 1998), globular clusters (Djorgovski & Meylan, 1994) and galactic nuclei 
(Schrodel et al. 2002). The first three columns list the cluster type, the total mass 
(in solar units) and the virial radius (in pc). For globulars the total mass and 
virial radius are given as distributions with a mean and the standard deviation 
around the mean. The orbital separation (in solar units) for a 1000 kT binary 
consisting of two 10 M stars is given in the fourth column. The fifth and sixth 
columns list the expected number of black-hole binaries that are formed by the 
cluster and the fraction of these binaries which merge within 12 Gyr (allowing 
~ 3 Gyr for the formation and ejection of the binaries and assuming a 15 Gyr old 
Universe (Freedman et al. 2001). The contributions to the total black hole merger 
rate per star clusters per year (MR) are given in the final column (for details 
The bottom row contains estimated parameters for the zero-age 
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population of globular clusters in the Galaxy. 

* Estimate for the parameters at birth for the population of globular clusters 
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Table 1.2. Galaxy morphology class, space densities, average absolute magni- 
tude (Heyl et al. 1997), and the specific frequency of globular clusters Sn (from 
van den Bergh, 1995 and McLaughlin 1999). The final column gives the contribu- 
tion to the total number density of globular clusters. The galaxy morphologies are 
identified as: E (elliptical), SO-Scd (spiral galaxies), Blue E (young blue elliptical 
galaxies), Sdm (blue spiral galaxies) and StarB (star burst galaxies). 



Galaxy 




</>GN 


M v 


S N h z 


GC space density 


Type 


[io- 


~ 3 ho Mpc" 3 ] 






[hi Mpc" 3 ] 


E-S0 




3.49 


-20.7 


10 


6.65 


Sab 




2.19 


-20.0 


7 


1.53 


Sbc 




2.80 


-19.4 


1 


0.16 


Scd 




3.01 


-19.2 


0.2 


0.03 


Blue E 




1.87 


-19.6 


14 


1.81 


Sdm/StarB 




0.50 


-19.0 


0.5 
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cluster M15. The globular cluster 47Tuc. is highly concentrated, but not 
quite in a state of core collapse. 






Figure 1.6. Three globular clusters to illustrate the effect of core collapse. 
Left: the rather shallow globular clusters Omega Cen. (NGC 5139, image by 
P. Seitzer) with a concentration Wo — 6. Middle: the concentrated cluster 
47 Tuc. (NGC 104, image by W. Keel) with concentration Wo ~ 9. Right: 
the core collapsed cluster M15 (NGC 7078, HST image R. Guhathakurta) with 
concentration Wo — 12. 



1.1.2 Fundamental time scales 

The evolution of a star cluster is dominated by two main effects; the mu- 
tual gravitational attraction between the stars and by the evolution of the 
individual stars. Before we discuss these two ingredients in more detail is it 
important do understand how they are coupled in the ecological network. 

1.1.2.1 Stellar evolution 

The fundamental timescale for stellar evolution is the nuclear burning 
timescale. This timescale only depends on characteristics of the stars them- 
selves and is unrelated to any of the dynamical cluster parameters 3 . The 
mass is the most fundamental parameter in the stellar evolution process; 
more massive stars burn-up more quickly than lower mass stars. The main 
sequence lifetime of a star with mass m and luminosity I can, to first order, 
be approximated with 

^LlxlO^ar^^f). (1.3) 

3 At least as long as the stars do not strongly influence each-other's evolution due to 
dynamical interactions. 
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After the main sequence the star generally grown to giant dimensions, 
but remains large only for approximately 10% of its main-sequence life- 
time . Massive stars continue their evolution after central hydrogen is 
exhausted by burning Carbon, Neon, Oxygen and Silicon until an iron core 
forms, which ultimately collapse catastrophically. The result is a supernova 
and the formation of a compact stars. Lower mass stars cannot process all 
material and shed their envelope in a planetary nebula phase. This results 
in a white dwarf. Rather detailed stellar evolutionary tracks are published 
in in the form of comprehensive look-up tables (Schaller et al 1993; Meynet 
et al. 1994) or fitting formulae (Eggleton, Fitchet & Tout 1989; Hurley et 
al 2000); using a variety of stellar evolution codes. 

1.1.2.2 Dynamical timescales 

The most fundamental dynamical timescale id for a star cluster is the 
crossing time or dynamical timescale, which for a cluster with half-mass 
radius R and dispersion velocity (v) can be written simply as 

t d = R/(v). (1.4) 

We can generalize this by using the cluster half-mass radius and the disper- 
sion velocity to calculate the global crossing time ih m of the star cluster. 

So long as stellar evolution remains relatively unimportant, the clus- 
ter's dynamical evolution is dominated by two-body relaxation, with pro- 
ceeds via the characteristic half-mass relaxation time scale (Spitzer, 1987) 



R*V /2 N 



GM 81nA 



(1.5) 



Here G is the gravitational constant, M, is the total mass of the cluster, 
N = M/(m) is the number of stars and r is the characteristic (half-mass) 
radius of the cluster. The Coulomb logarithm In A ~ ln(O.liV) = O(10) 
typically. In convenient units the two-body relaxation time becomes 

R \ 3/ V M \ 1/2 /lM m \ „ .,_! 



^""nisJ vm) lwT lnA) ■ (L6) 

Table fOl summarizes the relevant fundamental characteristics of var- 
ious types of star clusters. The parameters we selected are the stellar evo- 
lution time scale, cluster mass, size and velocity dispersion. With these we 
can compute the relaxation time and crossing time for the various clusters. 
Near the bottom of the table we present the collision time: t co \\ = l/r co n. 

4 As a rule of thumb, one can adopt that each subsequent burning stage for a star takes 
about 10% of the previous burning stage; Carbon burning lasts for about 10% of the 
helium burning stage, etc. 
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Here r co n is the collision rate for a cluster with stellar density n, and can be 
written as r co ii = nav, where a is the cross section for physical collisions: 



= ird 2 



v 
1 + — 

V 



(1.7) 



Here v^ is the relative velocity of two stars with mass mi and m-i at infinity 
and v is the relative velocity at closest distance d in a parabolic encounter, 
i.e: v 2 = 2G(m 1 + m 2 )/d. The second term results from the gravitational 
attraction between the two stars, and is referred to as gravitational focus- 
ing. 

We can now estimate the timescale for a collision between two stars as 
tco\\ = l/r co ii, which for a cluster with low velocity dispersions (woo <C v, 
can be written as 



t coll = 7 x 10 iu yr 



10 5 pc~ 3 \ / v, 



VlOkm/sy' ^ d 



M 



O 



for u»« 



(1-8) 

For a collision between two single stars we can adopt d ~ 3r* (Davies, Benz 
& Hills 1991). 

The last column in table ll.31 gives, in retrospect to the star clusters, a 
comparison with the beautiful city in which this meeting is organized. 



1.1.3 The effect of two-body relaxation: dynamical friction 

A star cluster in orbit around the Galactic center is subject to dynamical 
friction, in much the same way as dynamical friction drives massive stars 
toward the cluster center. This causes clusters to spiral into the Galactic 
center and stars to the cluster center. A star cluster is generally destroyed 
by the tidal field when it approaches the galactic center (see Gerhard 2001). 
We derive here the dynamical friction time scale for a mass point in the 
potential of the Galactic center. The derivation of the dynamical friction 
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time scale of a star as it spirals to the cluster center is very similar, with the 
major exception than the cluster potential is much more complicated that 
the potential of the Galaxy. We therefore opted for showing the derivation 
for a galaxy instead. 

We assume the inpiraling object to have constant mass M, deferring 
the more realistic case of a time-dependent mass (for example in the case 
of a star cluster sinking to the Galactic center, see Ea II. 3711 to McMillan 
& Portegies Zwart (2003). 

The drag acceleration due to dynamical friction in an infinite ho- 
mogeneous medium with isotropic velocity distribution that is not self- 
gravitating is (equation [7-18] in Binney & Tremaine, 1987) 



4vrlnAG 2 Af / 9 G (i? gc ) 



erf (X) - ^e~ x ' 

\/7T 



(1.9) 



Here In A is the Coulomb logarithm for the Galactic central region, for 
which we adopt In A ~ R gc /R, erf is the error function and X = Uc/v^fdisp, 
where i>disp is the one-dimensional velocity dispersion of the stars at dis- 
tance i?g C from the Galactic center, and v c is the circular speed of the 
cluster around the Galactic center. 

The mass of the Galaxy lying within the cluster's orbit at distance R gc 
( £, 500 pc) from the Galactic center is (Sanders & Lowinger 1972; Mezger 
et al. 1996) 

M G ai(i? g c) - 4.25 x 10 6 (^\ M . (1.10) 

Its derivative, the local Galactic density (see Portegies Zwart et al. 2001a) 
is 

PG (R gc ) ~ 4.06 x 10 5 (jZ£\ M Q pc- 3 . (1.11) 

For inspiral through a sequence of nearly circular orbits, the function 
erf (AT) — ^£ exp(— X 2 ) appearing in Eq. 11.91 may be determined as fol- 
lows. 

Following Binney & Tremaine (p. 226), we write the equation of dy- 
namical equilibrium for stars near the Galactic center as 

dP _ GAfdCiZge) 

where P = kTp/(m), \kT = \{m){v 2 ). Since u 2 = \{v 2 ), it follows that 
P = u 2 p, and Eq. II. 121 becomes 

Uu 2 p) = - P -v c 2 , (1.13) 

dr r 
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where v c is the circular orbital velocity at radius R: fdisp" = GMQ a \(R gc )/R gc . 
For Mcai oc R gc x (see Eg . 11.101 . and assuming that «disp 2 oc v c 2 ~ i? 
we find ■y d i sp 2 p < 



R gc 2x ~ 4 , so 



•AC 



r—(v d isp 2 p) = (2a: - 4)u disp 2 p = -pt> c , 
ar 

and hence X = y/2 — x. Eq. 11.91 then becomes 



(1.14) 
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IX 
erf (X) exp(-A 2 



(1.15) 



For 



1.2, X = 0.89 and 



a = -0.411nA 



GM 



R, 



2 ' 



(1.16) 



Again following Binney & Tremaine, defining L — R gc v c and setting 
dL/dt = aR gc , we can integrate Eg . 11.161 with respect to time to find an 
inspiral time from initial radius Ri of 



T{ 



1.28M Gal (fli) 



In A 



AI 



GM GaX (Ri 

R2 



-1/2 



K^rra*' 



(1.17) 
(1.18) 



For definiteness, we have assumed In A ~ 4 (A ~ R gc /R ~ 100) in Eq. 11.181 
corresponding to a distance of about 10-30 pc from the Galactic center. 

1.1.4 Simulating star clusters 

Stars move around due to their mutual gravity. This principle was first 
accurately described by Sir. Isaac Newton in 1687 in his Philosophiae nat- 
uralis principia mathematica (an excellent short biography can be found at 
http : //www-gap . dcs . st-and . ac . uk. 

Newton's equation describe the gravitational interaction between two 
stars with masses mi and m 2 with relative positional vector r = r 2 — r^, 
which we can written as: 



dr^ 
di 2 " 



mi +m 2 
Ct r. 



(1.19) 



The minus sign in the right-hand side indicates that the interaction force 
(dr 2 /dt 2 ) is attractive. 

This second order differential equation can be integrated in many ways. 
At this moment Hut and Makino are in the process of writing a series 
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of 10 books about integrating this equation using the so called direct N- 
body technique. The first three volumes of this series are available at 
http : //www . ArtCompSci . org. Other recent excellent work is published by 
Heggie & Hut (2003) and Aarseth (2003). Therefore instead of worrying 
about the intricacies of TV-body techniques we continue directly with the 
core problem. 

First, however, it may be useful to explain a bit about the methods 
under consideration and some of its alternatives; it is not my intention to 
give a thorough overview of all the way in which you can solve an iV-body 
system, but it is good to have some, overview. In a direct A^-body solver you 
integrate the equations of motion of all stars in the system by computing 
the forces from each star directly. This means that the amount of work for 
the computer scales roughly with the square of the number of particles N, 
or in units of CPU time: 

tcpu = (T) . (1-20) 



which for large N becomes tcpu = hmjv^ooiV 2 . This scaling becomes 
even worse if one imagines that the dynamical time unit in a star cluster 
is inversely proportional to N (see Eq. ll.4p . resulting in a time complexity 
which approaches icpu oc N 3 . Due to this high computational cost it 
is at this moment not possible to integrate the equations of motion of 
10 6 stars for a Hubble time. The largest simulations so far have been 
done are N ~ 10 5 stars for a Hubble time (Baumgardt et al 2003) and of 
N = 585.000 stars for the first 12Myr of the cluster (Portegies Zwart et al 
2004). 

One way to escape this conundrum is by integrating the equations of 
motion less accurately, or to not integrate them at all but by approximating 
the time evolution of the (grand)canonical ensemble of stars which forms 
the cluster. The (semi) approximate methods are generally harder to code 
and often impossible to assess, which makes fine-tuning to direct iV-body 
simulations inevitable. The core problem here is that a star cluster is in 
a delicate way not in perfect virial equilibrium. Methods which assume a 
vitalized state therefor have a distinct disadvantage over methods which 
do not explicitly require equilibrium. Often the more complex numerical 
coding of approximate solvers is well spend for large systems, like galaxies 
or cosmological simulations, in which low precision methods are preferred 
because of the shear number of 'stars' which have to be modeled. 

For simulating dense star clusters the best way is probably still the 
direct integration of the TV-body system, though competitive approaches 
have been taken (see e.g., § 11. 1.4. II and § 11.1. 4. 2JI . 

An extensive comparison between various types of N-body codes has 
been performed in Heggie's (1998) collaborative experiment, the results of 
which can be inspected at 

http: //www. maths . ed. ac . uk/^douglas/experiment .html. 
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Recently Spinnato et al. (2003) carry out a comprehensive comparison 
between three very different N-body techniques. The methods they adopt 
are a direct integration approach, which is, though accurate, strongly lim- 
ited in the number of particles which can be integrated. For larger particle 
numbers they used a tree-code, and for the same system but with up to 
several million particles they adopted a particle-mesh technique. 

1.1.4-1 Particle mesh 
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Figure 1.7. The different grids of SUPERBOX (Fellhauer et al. 2000) for 4 
cells per dimension. The finest and intermediate grids are focused on the object 
of interest. Figure from Spinnato et al (2003). 

To perform calculations for close-to homogeneous particle distributions 
a particle-mesh code is quite suitable. Grainier systems like dense star 
clusters are less suited, though some advances have been made by increasing 
resolution in substructure regions. However, caution has to be taken to 
make sure that the studied stellar system does not relax, as relaxation is 
generally not treated correctly in a particle-mesh technique. A nice example 
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is SUPERBOX (Fellhauer et al. 2000) in which accuracy is sacrificed for 
speed. 

In the particle-mesh technique densities are derived on Cartesian grids. 
Using a fast Fourier transform algorithm these densities are converted into 
a grid-based potential. Forces acting on the particles are calculated us- 
ing these grid-based potentials, making the code nearly collision-less. To 
achieve high resolution at the places of interest several techniques to im- 
prove a better local accuracy are used; SUPERBOX for example, incorpo- 
rates two levels of sub-grids which stay focused on the objects of interest 
while they are moving through the simulated area (see Fig. ll.7f) . providing 
higher resolution where required. 

Particle mesh codes, however, will always suffer from discrete effect 
due to the projection of the system on a grid. The size of a grid cell how- 
ever, appears to be related directly to the softening parameter e used in 
tree codes and in some direct A^-body codes. The concept of softening as 
introduced by (Aarseth 1963) is a technique to prevent large angle gravi- 
tational scattering by increasing the distance between two particles with a 
small parameter e. As long asr > e the effect of softening is not notice- 
able, but at close distances it has a profound effect on the behavior of the 
system. 

In order to understand how the cell length of the particle-mesh code 
and the softening parameter e of direct Af-body and treecodes relate with 
each other, we compare in Fig. II. 81 the results from the particle-mesh code 
SUPERBOX with the GADGET (Spingel et al 2001) treecode simulations 
for 80 000 particles. In this example Spinnato et al. (2003) use a black 
hole of mass 0.00053 (the mass of the inner part of the galaxy is unity in 
these units) to sink to the center of the Galaxy from a normalized distance. 
We can scale these numbers to astrophysically relevant units, in which case 
the black hole is ~ 65000 M and born at a distance of ~ 8 pc from the 
Galactic center. The initial orbit of the black hole was circular, but it still 
sinks slowly to the center of the Galaxy, due to dynamical friction (see 
§ 11.1.3(1 . Figur eTOl shows the distance of the black hole to the Galactic 
center as a function of time for the two computer codes. The results are 
presented for two values of the softening parameter e in the treecode and 
compared with two values of the cell-size in the particle mesh code. 

The in-fall of the black hole, as shown in Fig. 11.81 depends on the 
values of I and e in a remarkably similar way; I and e seem to play the 
same qualitative role, but also quantitatively the results are quite similar. 

1.1.4-2 Tree codes 

The concept of a hierarchical treecode was introduced by Apple (1985) and 
Barnes & Hut (1986), and is now widely used for the simulation of (near) 
collision-less systems. Figurc lPJl illustrates the concept of hierarchically 
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Figure 1.8. Spiral in of a black hole (of mass Mbh — 0.0005) to the center 
of the Galaxy (with mass MGai = 1)- The simulations ware performed with 
80,000 stars, using the same initial realization for both the treecode and the 
particle-mesh code. The particle-mesh simulations are run with cell size ~ 0.037 
(left) and 0.086 (right); softening parameters in the treecode runs are resp. 0.030 
(left) and 0.060 (right). 



deviding the spacial coordinates in the tree code. The force on a given par- 
ticle in a treecode is computed by considering particle groups of ever larger 
size as their distance from the particle of interest increases. Force contri- 
butions from such groups are evaluated by truncated multipole expansions. 
The grouping is based on a hierarchical tree data structure, which is real- 
ized by inserting the particles one by one into initially empty simulation 
cubes. Each time two particles are in the same cube, it is split into eight 
'child' cubes, whose linear size is one half of its parent's. This procedure is 
repeated until each particle is in a different cube. Hierarchically connecting 
such cubic cells according to their parental relation leads to the hierarchical 
tree data structure. The force on the particle of interest is then computed 
on neighboring cubes, which increase in size as they are further away. One 
of the interesting characteristics of tree-codes is the relatively simple par- 
allelization by domain decomposition of the spacial coordinates (Olson & 
Dorband 1994); a disadvantage is the lack of support for special purpose 
hardware (see however Kawai et al 2004). 

In the same way we compared a particle-mesh code with a treecode 
in § 11.1.4.11 we can also compare the tree code with a direct TV-body code 
(see next section). In Fig. ll.lOl such comparison is illustrated on the same 
simulation of the black hole which spirals to the Galactic center. Both 
codes are run with N = 80 000 and for a softening of e ~ 0.0037, 0.030 and 
e ~ 0.060. In this comparison we adopted kira from the Starlab package 
as the direct TV-body code. The latter code was also run with zero softening 
(e = 0) but for the treecode this did not produce reliable results. 

Interestingly, from the comparison between the particle-mesh code, 
the tree code and the direct TV-body code in the figures lTTHl and n~TUI it is 
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Figure 1.9. Schematic illustration of tree-building in a Barnes & Hut (1986) 
algorithm in two dimensions (from Springel, Yoshida and White, 2001). The 
particles are first enclosed in a single square, which is iteratively subdivided in 
squares of only half the size until a single particle remains per square. 
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Figure 1.10. Comparison of results from the PP code with results from the 
treecode, at different values of e. For all cases shown here is N — 80 000 and 
Mbh — 0.000528. The PP simulation with e = 8eo has been already shown in 
Fig. ??. the runs ware performed with e ~ 0.0037 (left), 0.030 (middle) and 
e ~ 0.060 (right). 
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evident that in this practical case all three codes produce qualitatively and 
quantitatively the same results. One can wonder why then it is so important 
to perform an accurate calculation, as low resolution simulations produce 
results which are consistent with high precision simulations? 5 

I.I.4.3 Direct N -body 

The most most accurate way to simulate the dynamical evolution of a star 
cluster is by solving Eg . 11. 191 using direct integration. Direct iV-body codes 
seem simple to write, as it is just a matter of solving Eg . 11. 191 in small 
steps, but in fact a computer program that does the job accurately and 
quickly is very hard to write right. This technique was pioneered by von 
Hoerner (1960, 1963) Aarseth (1963), Aarseth & Hoyle (1964) and van 
Albada 1968). Excellent reviews are published in the earlier mentioned 
books by Hut & Heggie (2003), Aarseth (2003) and Hut & Makino (2003), 
but see also Aarseth & Lecar (1975). A further reference to Heggie & 
Mathieu (1986) cannot be omitted as this valuable paper discusses the 
dimension less units used in A^-body techniques, in which M = R = G = \. 

In a direct N-body code the forces between all stars are calculated with 
numerical precision. The time complexity in this calculation is, as discussed 
earlier 0(iV 3 ). In such a code, the particle motion is followed using a 
high-order integrator often with a predictor-corrector scheme (Makino and 
Aarseth 1992). These codes can generally not work with shared time steps 6 ; 
it saves time and gains accuracy to allow stars in a strong encounter to be 
updated by individual time steps 7 . For simpler parallelization one generally 
adopts block time steps so that groups strongly interacting particles are 
integrated more frequently than weakly interacting stars (McMillan 1986a; 
1986b; Makino 1991). Still special treatment for binaries and higher order 
hierarchies are required to prevent the code to come to a grinding halt 
during strong encounters. 

During a time step, particle positions and velocities are first predicted 
to fourth order using the acceleration and "jerk" (time derivative of the 
acceleration), which are known from the previous step. The new acceler- 
ation and jerk are then computed, and the motion is corrected using the 
additional derivative information. 

One of the great advantages of using a direct TV-body solver is the 
simplicity at which extra effects can be incorporated. Since each star in the 
cluster is represented by a particle in the code, individual characteristics, 
such as stellar properties, can be accounted for relatively easily and without 

5 It is probably worth revisiting this problem by performing a details side-by-side anal- 
ysis of two or more iV-body simulation methods (see Heggie et al 1998 and his collabo- 
rative experiment at http://www.maths.ed.ac.uk/~douglas/experiment.html). 

6 All particles share the same time step. 

7 Each particle is integrated with its own special timestep, particles in a strong inter- 
action can then be integrated more accurately than weakly interacting particles. 
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loss of generality. This makes the direct TV-body method preferable for 
simulating star clusters where these effects are important. In our case 
we are interested in the evolution of black holes, which are relatively rare 
objects. It would therefore be best to utilize a technique in which we can 
also treat the black holes individually. The draw back here is that black 
holes are so rare that large clusters have to be simulated in order to obtain 
enough statistics on the black hole population. 

On the other hand the gravitational N-body problem has many appli- 
cations over a wide range of research fields, including informatics, compu- 
tational science, geology and astronomy. 

1.1.4-4 GRAPE family of computers 

The enormous computational requirement for solving the ./V-body problem 
with the direct method has been effectively addressed by a small team of re- 
searchers, who developed the GRAPE family of special purpose computers. 
GRAPE (short for GRAvity PipE) hardware was designed and built by a 
group of astrophysicists at the University of Tokyo (to name only the most 
relevant publications: Sugimoto et al. 1990; Fukushige et al 1991; Ito et al. 
1991; Okumura et al 1993; Taiji et al 1996; Makino et al 1997; Kawai et al. 
2000; Makino 2000; Makino et al. 2003). It may be clear that GRAPE is 
a very successful endeavor. The history and computational science of the 
GRAPE project is published by Makino & Taiji (1998). 

The GRAPE family of computer are like a graphics accelerator speed- 
ing up graphics calculations on a workstation, without changing the soft- 
ware running on that workstation, the GRAPE acts as a Newtonian force 
accelerator, in the form of an attached piece of hardware. In a large- 
scale gravitational N-body calculation, where N is the number of particles, 
almost all instructions of the corresponding computer program are thus 
performed on a standard workstation, while only the gravitational force 
calculations, the innermost loop, are replaced by a function call to the 
special-purpose hardware. 

Figur enTTl shows the fully configured GRAPE-6 at Tokyo University. 

1.1.5 Performing a simulations 

Before starting to simulate one may want to consider what technique is 
most suitable. In our further discussion that will be the direct iV-body 
integrator. Several such computer programs are readily available. The 
starlab environment provides a entire library of functions and routines 
built around the main ./V-body integrator. The package can be downloaded 
from 

http : //www/manybody . org/starlab . html. But also NBODY 1-6 are avail- 
able by ftp via ftp : //ftp . ast . cam . ac . uk/sverre/. Both codes are large 
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Figure 1.11. The large 64 Tflops GRAPE-6 configuration in Tokyo, in the 
summer of 2003. 



and very complicated as they have been evolving to the current sophistica- 
tion over more than a decade. But simpler alternatives are available from a 
variety of sources (from example via: 

http : //www . ids . ias . edu/^piet/act/comp/algorithms/starter/. A fully 
operational parallel A^-body code based on the above mentioned starter 
code can be obtained from 

http: //carol . science . uva.nl/^spz/act/modesta/Sof tware/index.html. 

Whatever computer code you select or even if you write one from 
scratch, make sure that you test it. Test the code against other similar 
codes, test it with calculations by hand, regardless how painstaking this 
often is, and test it against simple problem for which the solution is known. 
You should develop a feeling of the regimes where the code can be trusted 
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and in what cases extra care must be taken in interpreting the results. 

Let's assume that we have found an interesting problem, which we 
think we can solve with the code available. The main problem of starting 
a simulation then is the selection of initial conditions. Starting with wrong 
initial conditions is a complete waste of time. It is better to spend enough 
time thinking about the initial conditions, until you are convinced that they 
are the best choice. Possibly you want to perform several test calculations 
to converge to a better understanding of the question asked and the initial 
conditions required to give the most reliable answer. 

For simulating a star cluster the primary initial conditions are: 

• What are the basic cluster properties: mass, size, density profile? 

• How many stars did the cluster have at birth? 

• How are the stellar masses distributed? 

• What is the fraction of primordial binaries? 

• And what are the binary parameters: semi-major axis, eccentricity, 
inclination, etc. 

• Do you want to include triples and higher order systems? 

• and what about the shape and strength of the external tidal field of 
the Galaxy? 

• Are there tidal shocks, spiral arms or other external potentials to worry 
about? 

• Is there anything else to add like: passing molecular clouds or black 
holes, etc. 

Many of the above effects have to some extend been incorporated in 
various calculations. And there is a rich scientific literature about the 
relevance and effect of many of these ingredients. 



1.2 Theory of star cluster evolution 

Now we have set the stage and discussed the tools and the techniques we 
can continue by discussing the global evolution of star clusters, which is 
characterized by three quite distinct phases; these are subsequently: A the 
early relaxation dominated phase, followed by phase B in which the ~ 1 % 
(by number) most massive stars quickly evolve and lose an appreciable 
fraction of their mass. Finally, phase C starts when stellar evolution slows 
down even more and relaxation takes over until the cluster dissolves, to 
complete the list we can define phase D which is associated with the final 
dissolution of the cluster due to tidal stripping, but we will not discuss this 
phase in detail. 
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1.2.1 Phase A: t < lOMyr 

In an early stage, when stellar evolution is not yet important the star cluster 
is dominated by its own dynamical evolution; we call this phase A. This 
stage in the evolution of the cluster is only relevant if i r t ^S 100 Myr. For 
most open star clusters and for all globular clusters this stage is probably 8 
not very important, but for YoDeCs it is crucial, as we explain below. 

In the following discussion we assume, for clarity, that the location in 
the cluster where the stars are born is unrelated to the stellar mass, i.e.: 
there is no primordial mass segregation. In that case, the early evolution of 
the star cluster is dominated by two-body relaxation, or to be more precise 
by dynamical friction. 

In young dense clusters dynamical friction implies a characteristic time 
scale tdf for a massive star in a roughly circular orbit to sink from the half- 
mass radius R, to the cluster center (Spitzer 1971; 1987): 

(m) 0.1387V /i? 3 \ 1/2 

df lOOM ln(O.HM/lOOM ) \GM ) ' V ' ' 

Here (m) and M are the mean stellar mass and the total mass of the cluster, 
respectively, N is the number of stars, and G is the gravitational constant. 
For defmiteness, we have evaluated tdf for a 100 M Q star. Less massive 
stars undergo weaker dynamical friction, and thus must start at smaller 
radii in order to reach the cluster center on a similar time scale. 

Dynamical friction will have two very distinct effects on the cluster: 
1) it tends to produce cores in hitherto core-less clusters, and 2) it initiates 
core collapse in other clusters. These two statements seem to contradict 
each other but as we will see below, this is not the case. 

1.2.1.1 Dynamical friction induced core development 

Consider a gravitationally bound stellar system in which most of the mass 
is in the form of stars of mass m ~ (m) , but which also contains a subpop- 
ulation of more massive objects with masses Wh- The orbits of the more 
massive objects decay due to dynamical friction. Assume that the stellar 
density profile is initially a power-law in radius, p(r) oc (r / a) - " 1 M / a with 
density scale length a (Dehnen 1993; Merritt et al. 1994). The orbits 
of the massive stars decay at a rate that can be computed by equating 
the torque from dynamical friction with the rate of change of the orbital 
angular momentum. We adopt the usual approximation (Spitzer 1987) in 
which the frictional force is produced by stars with velocities less than the 
orbital velocity of the massive object. The rate at which the orbit decays, 

8 Regretfully we know very little about the early stages of globular clusters, and it is 
therefore hard to say whether phase A was important. 
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assuming a fixed and isotropic stellar background, is 



§~^V™S-© 7/M ™ <-> 



with 
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2-7, 

(1.23) 

Here (3 = (6 — 7)/2(2 — 7) and In A is the Coulomb logarithm, roughly equal 
to 6.6 (Spinnato et al. 2003). For 7 = 1.0 (2.0), F = 0.19 (0.43). 

If we approximate the cluster structure with an isothermal sphere, 
we find (Binney & Tremaine 1987, Eq. 7-25) that a star of mass mj, at 
distance R from the cluster center drifts inward at a rate given by 

R dR On* 

dt (v) K ' 

Here (v) is the clusters' velocity dispersion. 

Eauation ll.221 implies that the massive object comes to rest at the 
center of the stellar system in a time 



M /i^ (6 -" )/2 



*« - 2 \/77T7 — ( L25 ) 

V GM m h \ a J v ; 

with Ri the initial orbital radius. 

Or, we can express the dynamical friction time in terms of the half- 
mass relaxation time by substituting Eg. II. 61 in Eg . 11. 241 and integrate with 
respect to time. 

tdf ~ 3.3^-trt . (1.26) 

m h 

To estimate the effect on the stellar density profile, consider the evo- 
lution of an ensemble of massive particles in a stellar system with initial 
density profile p ~ R~ 2 . The energy released as one particle spirals in from 
radius Ri to Rf is 2m\ 1 a 2 \n(Ri/Rf), with a the ID stellar velocity dis- 
persion. Decay will halt when the massive particles form a self-gravitating 
system of radius ~ GM^/u with Mjj = ^2 TO h- Equating the energy re- 
leased during in-fall with the energy of the stellar matter initially within 
r c , the "core radius," gives 

r ^^^ 1r {gm-J- (l27) 

Most of the massive particles that deposit their energy within R c will come 
from radii R4 sw a few x R c , implying R c « several x GM^/a 2 and a 
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displaced stellar mass of ~ several x M^ (see also Watters 2000). If M^ « 
10~ 2 M (Portegies Zwart & McMillan, 2000) then R c /a « several x 2M h /M 
and the core radius is roughly 10% of the effective radius. Merritt et al 
(2004) discuss this process in more detail and apply it successfully to the 
evolution of the core radii of large Magellanic cloud star clusters. 

Evolution will continue as the massive particles form binaries and begin 
to engage in three-body interactions with other massive particles. These 
superelastic encounters will eventually lead to the ejection of most or all 
of the massive particles. Assume that this ejection occurs via many small 
'kicks', such that almost all of the binding energy so released can find its 
way into the stellar system as the particle sinks back into the core after 
each ejection. The energy released by a single binary in shrinking to a 
separation such that its orbital velocity equals the escape velocity from the 
core is ~ m^a 2 ln(4Mh/M) (see also 8 I1.3.4JI . If all of the massive particles 
find themselves in such binaries before their final ejection and if most of 
their energy is deposited near the center of the stellar system, the additional 
core mass will be 

M c wM h ln(— J (1.28) 



M h 

e.g. ~ 5-Mh for M^/M = 0.01, similar to the mass displaced by the initial 
in-fall. The additional mass displacement takes place over a much longer 
time scale however and additional processes (e.g. core collapse) may com- 
pete with it. 

1.2.1.2 Dynamical friction induced core collapse 

The dynamical evolution of the star cluster drives it toward core collapse 
(Antonov 1962; Spitzer & Hart, 1971a; 1971b) in which the central density 
runs away to a formally infinite value in a finite time. In an isolated 
cluster in which all stars have the same mass, core collapse occurs in a 
time t cc ~ 15t rt (Cohn 1980; Makino 1996; Joshi et al 2001). 

If the dynamical friction time scale of a star cluster is shorter than the 
lifetime of the most massive stars, the cluster may experience an early phase 
of core collapse before the first supernova occurs. This can only happen if 
the initial half-mass relaxation time of the cluster is small t It ^ 100 Myr 
otherwise the most massive stars burn-up before they reach the cluster 
center (Portegies Zwart et al 1999; Giirkan et al 2004). The early core 
collapse in a cluster with small relaxation time is illustrated in figure HTT^l 
where we plot the core radius as a function of time for a number of young 
star clusters with a relaxation time of about 100 Myr. 

The details of what exactly happens in the cluster core, and whether 
the cluster will experience gravothermal oscillations probably depends quite 
critically on the initial density profile, as we will discuss in more detail in 
Scct. lDOl 
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Figure 1.12. Evolution of the core radius, defined as the radius at which the 
projected density falls to one-half of its central value. Each curve is the average 
of the various runs at the specified value of 7, 7 = 1 for the top curve, 1.5 for the 
middle curve and 7 = 2 for the bottom curve (see 5 11.2, 1.1^ . The figure is taken 
from Merritt et al. 2004). 



When after core collapse the most massive stars explode the core radius 
remains highly variable, but small on average. What exactly happens at 
this stage is still ill understood. Naively one would expect the collapsed 
core to expand, as stellar mass loss drives an adiabatic expansion of the 
cluster. This can be been seen in the post core collapse evolution of the 
bottom solid curve in figurc tTT^I Stage B starts when dynamical friction 
and relaxation cannot further drive the core collapse of the cluster, but 
expansion by stellar mass loss starts to dominate. Typically this happens 
around ~ 10 Myr. 

The cluster simulated for Fig. ll.l^l was MGG11, in the starburst galaxy 
M82. It contained 131072 single stars from a Salpeter initial mass function 
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between 1M Q and 100 M distributed in a King (1966) model density 
profile with dimension less dept Wo = 3 (shallow) to Wo = 12 (very con- 
centrated). The half mass radius of these simulated clusters was 1.2 pc. 

We discuss each curve in figure [TTT31 in turn, starting at the top. The 
core radius of the shallow model, c ~ 0.67 (Wo = 3), hardly changes with 
time. The intermediate model c ~ 1.8 (Wo = 8) almost experiences core 
collapse near t = 3 Myr but as stellar mass loss starts to drive the expansion 
of the core it never really experiences collapse. This is the moment where 
phase B sets in. Core collapse occurs in the c ~ 2.1 (Wo = 9) model 
near t — 0.8 Myr. The c ~ 2.7 (Wo = 12) simulation is so concentrated 
that it starts virtually in core collapse, and the entire cluster evolution is 
dominated by a post-collapse phase. At this point it is not a-priory clear 
why core radius for models Wo = 9 and Wo = 12 fluctuate so much more 
violently than in the models with smaller concentration. 

1.2.2 Phase B: lOMyr > t < lOOMyr 

After the first few million years and until the most massive stars have 
turned into compact remnants, the cluster will be dominated by stellar 
mass loss. Since the most massive stars evolve first and sink most quickly 
to the cluster center, mass is lost from deep inside the potential well. A 
collapsed cluster may recover from its earlier core collapse due to stellar 
mass loss (see S ll.2.1.2|) . The dotted curve in fig. II. 131 (for the simulations 
with Wo = 12) illustrates the growth of the core as a result of heating by 
dynamical friction and stellar evolution mass loss is quite clearly visible; it 
is hard to separate the two effects as both are taken into account in the 
calculation self consistently. For the slightly shallower initial density profile 
(Wo = 9) the steady growth of the core radius after collapse is less clear, 
but the rate is similar as for Wo = 12. This indicates that it is indeed 
mainly stellar mass loss which drives the expansion. At this stage we do 
not know why the core radius in the Wo = 9 models fluctuate more wildly 
than in the Wo = 12 simulations. Details probably depend quite sensitively 
on the presence of an intermediate mass black hole, which could form in the 
proceeding phase A. Such massive compact object can effectively heat the 
cluster core as it forms tight multiple systems with other (massive) stars 
(Baumgardt et al 2003). 

Few studies has been carried out for clusters with short relaxation time 
to understand this particular evolutionary stage. For a long relaxation time 
i r t ~ 100 Myr, it is quite clear that the cluster expands substantially during 
the first Gyr of its lifetime (Chernoff & Weinberg 1990; Fukushige & Heggie 
1995; Takahashi & Portegies Zwart 2000; Baumgardt & Makino 2003). 

The evolution of the cluster mass if given in Fig. II. 14l ( taken from Taka- 
hashi & Portegies Zwart 2000) for a variety of particle numbers, ranging 
from 1024 to 32768. The first few 100 Myr of all clusters is very similar, but 
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Figure 1.13. Evolution of the core radius for four simulations of the star 
cluster MGG-11. These calculations are performed with Starlab with c ~ 0.67 
(Wo = 3), c ~ 1.8 (Wo = 8), c ~ 2.1 (Wo = 9, bold curve) and c ~ 2.7 (Wo = 12, 
dotted curve) indicated along the right edge of the figure as W3, W8, W9, and 
W12 respectively. The W9 curve is plotted with a heavy line to distinguish it 
from the curves for W8 and W12. The age range of the observed clusters MGG11 
is indicated near the bottom of the figure. In figur eH^TI we present the evolution 
of the average core density for the simulations with Wo = 3, 8 and 9. 



at later time large differences appear. In the first epoch, during the first 
few hundred million years, about 20% of the cluster mass is lost. This mass 
loss is the result of stellar evolution and, in less extend by tidal stripping. 
Tidal stripping and relaxation become important at later time. 



1.2.3 Phase C: t > lOOMyr 

After the most massive stars have turned into remnants, stellar evolution 
slows down, and relaxation processes can take over again. In fie. II. 141 this 
phase starts at an age of about a Gyr. The reason that stellar mass loss 
slows down is two-fold, 1) low mass stars remain on the main sequence 
longer than high mass stars, and 2) once the star turns into a remnant, 
the mass lost in the process is relatively small; imagine a 12 M Q star loses 
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Figure 1.14. Mass as a function of time for a number of Fokker-Planck models. 
The four solid lines represent the results of model Aa with 32k, 16k, 4k and Ik 
particles from left to right, respectively. Dotted curves present model Ie for the 
same numbers of particles as for model Aa (model names refer to Takahashi & 
Portegies Zwart, 2000). The difference between the models Ie and Aa hides in 
the way escaping stars are identified. 



about 10.6 Mq by stellar wind and in the supernova explosion, whereas a 
2M turns into a 0.64 Mq white dwarf, losing only 1.36 M q in the process. 

The effects of mass lost from the evolving stars and mass lost from the 
dynamical evolution of the cluster are coupled. If stellar mass loss slows 
down, the cluster responses to this by a slower expansion, which again 
makes it less prone to tidal stripping. While the importance of stellar 
evolution diminishes, relaxation gradually takes over until it becomes the 
dominant mechanism which drives the evolution of the cluster. 

The later stage of the low N (Ik and 2k) clusters in Fig. 11.141 are much 
stronger affected by relaxation than the high N (16k and 32k) clusters. This 
effect was named the ski-jump problem in Portegies Zwart et al. (1998). 
The transition between ski-clusters (low N in fig. ll.14p and jump-clusters 
(high N) is a result of the non-linear interaction between the external tidal 
field of the parent galaxy and relaxation. 

At later time during phase C dynamical friction once again become 
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important. This time not driven by massive stars, as these have all gone 
supernova by now, but by the compact remnants formed in supernovae; 
black holes and neutron stars, but also heavy white dwarfs, blue stragglers 
and giants. All these stellar species are generally more massive than the 
mean mass, and therefore subject to dynamical friction. A similar process 
as in phase A starts again and the cluster may experience core collapse for 
the second time. It may therefore be possible that a cluster experiences 
two very distinct phases of core collapse, one during phase A, and again at 
a later time, during phase C (see also Deiters & Spurzem 2000, 2001). 

In a realistic cluster, however, there are a number of additional com- 
plications which are particularly important at this stage, in part because 
it may take rather long to reach a state of core collapse again because the 
stellar mass function is rather flat now, with a relatively small difference 
between the least and the most massive stars. The consequence is that 
external influences, like disc shocks, passing molecular clouds and the pres- 
ence of an external tidal field, may become particularly important at this 
stage, simply because they have a lot of time to accumulate their effect. 

The slow-down of stellar evolution has a second important conse- 
quence, which is the termination of active binary evolution. Only relatively 
low mass stars are able to evolve off the main sequence, and no supernova 
will occur after ~ 10 8 years. It becomes therefore almost impossible to 
ionize hard binaries, which may effectively arrest the collapse of the clus- 
ter core (Fregeau et al. 2003). The binaries may therefore once more 9 
become dynamically important, and heat the cluster by interacting with 
single stars or other binaries (Heggie 1975). 

We illustrate phase C in Fig. ll,15l with a simulation of 10,000 identical 
point masses initially distributed in a Plummer sphere. Fig. ll.15l shows 
the evolution of the core radius of this model. Core collapse occurs at 
t cc ~ 15.2±0.1 £ r t (for these initial conditions t T t — 150 iV-body time units). 
This result is consistent with earlier calculations of e.g., Cohn (1980) and 
Makino (1996). Doubling the mass of 20% of the stars reduces the core 
collapse time to t cc ~ 7.2 i rt . Making 20% of the stars 10 or 100 times 
more massive reduced the time of core collapse further, to t cc ~ 1.4 t rt and 
t cc ~ 0.16 trt 5 respectively. Clearly, the more massive stars drive the core 
collapse of the cluster by dynamical friction, as idf — ^^■t c t (see Ea . 11.211 
and also Watters et al. 2000). 

The presence of a population of relatively massive stars, such as black 
holes will therefore shorten the timescale for core collapse in phase C of 
the cluster evolution, but even in the absence of black holes or other heavy 
remnants core collapse cannot be prevented. In phase A this role was 
played by the most massive main-sequence stars. 

Note that if it is the population of dark remnants, black holes, neutron 

9 The first time binary interactions ware relevant was during phase A. 
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Figure 1.15. Evolution of the core radius for two star clusters well beyond the 
moment of core collapse. Time (in iV-body units) is along the x-axis and radius 
(log r, also in TV-body) units along the y-axis. Initially each cluster contained 
10,000 single stars distributed in a Plummer sphere with the same realization of 
positions and velocities. The solid curve was computed with equal masses for all 
stars, the dashed curve is computed with 20% of the stars being twice as massive, 
these calculations ware performed without stellar evolution and the units along 
the axis are in dimension less -/V-body units. 



star and white dwarfs, which experience the core collapse, the lighter stars 
will not necessarily follow. So, the cluster may physically be in a state of 
core collapse, where the optical observer would measure a King profile (see 
Baumgardt et al 2003a; 2003b) 



1.2.4 The consistent picture 

In this section we have seen that a star cluster experiences three very dis- 
tinct evolutionary phases, each of which is dominated either by relaxation 
or by stellar mass loss. 

Figure Tl . 1 61 summarized these phases in a single simulation which con- 
tains all relevant physics. The simulation presented here is carried out to 
illustrate the formation mechanism for intermediate mass black holes in 
dense young star clusters. It was performed with 131072 stars from a 
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Salpeter initial mass function, 10% (13107) of the stars form a hard binary 
system with a close companion. The initial conditions for this simulation 
are explained in more detail in § 11.3.31 Here we only use the result of the 
simulation to illustrate the three distinct evolutionary phases in star cluster 
evolution, as with the adopted parameters each phase is clearly present. 

The three phases, A, B and C, are identified with the horizontal bars in 
fig. II. 161 It it not a-priori clear when one phase stops and the next starts, 
and some gray area has to be allowed in which both, stellar evolution 
and relaxation, may temporarily have similar effect on the cluster. These 
runs ware continued till 100 Myr. 10 Since this simulation was performed in 
isolation, the dissolution of the cluster will take quite a while (Baumgardt, 
Hut & Heggie 2002) 

1.3 Black holes in star clusters 

In the previous sections we have set the stage for the evolution of star 
clusters, and we have introduced the fundamental physics. We will now 
continue with the evolution of black holes and their progenitors in star 
clusters. 

1.3.1 The formation of intermediate mass black holes in phase 
A clusters 

Young star clusters, with a half mass relaxation time t rt ^ 100 Myr are, as 
we discussed in Sect. ll.2"Tl prone to dynamical friction, and therefore are 
likely to experience core collapse before stellar mass loss drives the expan- 
sion of the cluster. Realistic clusters have a broad range in initial stellar 
masses, generally from m min £ 0.1 M to m max £ 100 M . Adopting such 
a mass function as a condition at birth, the mean mass (to) ranges then 
from (to) ~ 0.39 M (Salpeter 1955) to about 0.65 M (Scalo 1986), de- 
pending on the specific mass function adopted. Here I like to stress that 
there is a large variety of initial mass function available, apart from the 
two above mentioned there are the Miller & Scalo (1979) , Koupa, Tout & 
Gilmore (1990) with some adjustments for high mass stars Kroupa & Wei- 
dner (2003, see also Kroupa, 2001). These mass functions seem to differ 
quite substantially, but for the dynamical evolution of the cluster they do 
not make a big difference, as for this process it is the ratio of most massive 
star to the mean mass which counts (see Ea. ll.21|l . 

In a multi-mass system, core collapse is driven by the accumulation of 
the most massive stars in the cluster center. This process takes place on 

10 There was no particular reason to terminate this calculation at about 100 Myr, other 
than that I got a bit tired of babysitting the run after three months straight on the 
GRAPE-6 at the University of Tokyo. In time, this cluster will experience core collapse 
again, to dissolve eventually. 
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Figure 1.16. Evolution of the core radius (lower line) 10% and the 25% La- 
grangian radii for a 128k star cluster (axis are in iV-body units). The various 
stages of cluster evolution are indicated near the bottom of the figure. In the 
early evolution of the cluster expansion of the core is driven by mass segregation, 
followed by phase B in which the cluster expands by stellar evolution. Near the 
end of the simulation the cluster tends to experience core collapse (Phase C). 



a dynamical friction time scale (Eg . 11.21)) . Empirically, wc find, for initial 
mass functions of interest here, that core collapse (actually, the appearance 
of the first persistent dynamically formed binary systems) occurs at about 
(Portegies Zwart & McMillan 2000; Fregeau et al 2002) 



tc. 



0.2 tr 



(1.29) 



This core collapse time is taken in the limit where stellar evolution is unim- 
portant, i.e. where stellar mass loss is negligible and the most massive stars 
survive until they reach the cluster center, i.e.: what we called phase A in 

§EZH 

Experimentally wc find that starting with a Scalo (1986) initial mass 
function and a King (1966) Wq = 6 density distribution the time to reach 
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the first core collapse is t cc ~ 0.19 ± 0.08 t Tt . Some of our simulations 
ware performed with high Wo ^ 9 concentration. Core collapse in these 
models occurred within about million years, which for these simulations 
corresponded to i cc ~ 0.01 £ r t, or about 7.5% of the core relaxation time. 
These findings are consistent with the Monte-Carlo simulations performed 
by Joshi, Nave & Rasio (2001). 

The collapse of the cluster core may initiate physical collisions between 
stars. The product of the first collision is likely to be among the most 
massive stars in the system, and to be in the core. This star is therefore 
likely to experience subsequent collisions, resulting in a collision runaway 
(see Portegies Zwart et al. 1999). The maximum mass that can be grown 
in a dense star cluster if all collisions involve the same star is M run , where 

— ;t^ =-A/"coii£racoii- (1.30) 

at 

Here A/" co ii and <5m co n are the average collision rate and the average mass 
increase per collision (assumed independent). We now discuss these quan- 
tities in more detail. Interestingly enough, Giirkan et al. (2004) performed 
comparable calculations with a Monte-Carlo TV-body code, which produces 
qualitatively the same results. Also the calculations of Portegies Zwart et 
al (2004) who used two independently developed TV-body codes (NB0DY4 
and Starlab), obtain similar results. 

1.3.1.1 The collision rate Af co ii 

A key result from the simulations of Portegies Zwart et al. (1999) was the 
fact that collisions between stars generally occur in dynamically formed 
("three-body") binaries. The collision rate is therefore closely related to 
the binary formation rate, which we can estimate from first principles. 

The flux of energy through the half-mass radius of a cluster during 
one half-mass relaxation time is on the order of 10% of the cluster poten- 
tial energy, largely independent of the total number of stars or the details 
of the cluster's internal structure (Goodman 1987). For a system without 
primordial binaries this flux is produced by heating due to dynamically 
formed binaries (Makino & Hut 1990). It is released partly in the form of 
scattering products which remain bound to the system, and partly in the 
form of potential energy removed from the system by escapers recoiling out 
of the cluster (Hut & Inagaki 1985). Makino & Hut argue, for an equal- 
mass system, that a binary generates an amount of energy on the order of 
10 2 fcT via binary-single-star scattering (where the total kinetic energy of 
the stellar system is |iVfcT). This quantity originates from the minimum 
binding energy of a binary that can eject itself following a strong encounter. 
Assuming that the large-scale energy flux in the cluster is ultimately pow- 
ered by binary heating in the core. It follows that the required formation 
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rate of binaries via three-body encounters is 

r bf ~icr 3 — . (i.3i) 

tit 

For systems containing significant numbers of primordial binaries, which 
segregate to the cluster core, equivalent energetic arguments (Goodman 
& Hut 1989) lead to a similar scaling for the net rate at which binary 
encounters occur in the core. 

The above arguments apply to star clusters comprising identical point- 
masses. In a cluster with a range of stellar masses, three-body binaries gen- 
erally form from stars which are more massive than average. After repeated 
exchange interactions, the binary will consist of two of the most massive 
stars in the cluster. Conservation of linear momentum during encounters 
with lower mass stars means that the binary receives a smaller recoil ve- 
locity, making it less likely to be ejected from the cluster. The binary must 
therefore be considerably harder — x^ 10 3 /cT — before it is ejected following 
a encounter with another star (see Portegies Zwart & McMillan 2000). 

However, taking the finite sizes of real stars into account, it is quite 
likely that such a hard binary experiences a collision rather than being 
ejected. A strong encounter between a single star and a hard binary gen- 
erally results in a resonant interaction. Three stars remain in resonance 
until at least one of them escapes, or a collision reduces the three-body sys- 
tem to a stable binary. For harder binaries it becomes increasingly likely 
that a collision occurs instead of ejection (McMillan 1986; Gualandris et 
al. 2004; Fregeau et al. 2004). In the calculations of Portegies Zwart et 
al. (1999) most binaries experience a collision at a binding energy of order 
10 2 fcT, considerably smaller than the binding energy required for ejection. 
Accordingly, we retain the above estimate of the binary formation rate (Eq. 
I1.31J) and conclude that the collision rate per half-mass relaxation time is 

Kda ~ I0~ 3 f c —. (1-32) 

Crt 

Here we introduce f c < 1, the effective fraction of dynamically formed 
binaries that produce a collision. Note again that Eg . 11. 321 is valid only in 
the limit where stellar evolution is unimportant. 

The most massive star in the cluster is typically a member of the inter- 
acting binary and therefore dominates the collision rate. Subsequent colli- 
sions cause the runaway to grow in mass, making it progressively less likely 
to escape from the cluster. The star which experiences the first collision is 
therefore likely to participate in subsequent collisions. The majority of col- 
lisions thus involve one particular object — the runaway merger — generally 
selected by its high initial mass and proximity to the cluster center. 

For systems containing many primordial binaries the above argument 
must be modified. Since dynamically formed binaries tend to be fairly 
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soft — a few kT — the fraction of interactions with primordial binaries lead- 
ing to collision is comparable to the value f c above. However, a critical dif- 
ference is that, in systems containing many binaries, the collisions involve 
many different pairs of stars, not just the binary containing the massive 
runaway. The total collision rate is therefore much higher, but most colli- 
sions do not contribute to the growth of the runaway merger. The presence 
of primordial binaries have little influence on the collision runaway. 11 

1.3.1.2 Average mass increase per collision 

Once begun, the collision runaway dominates the collision cross section. 
The average mass increase per collision depends on the characteristics of 
the mass function in the cluster core. A lower limit for stars which par- 
ticipate in collisions can be derived from the degree of segregation in the 
cluster. Inverting Ea. ll.2ll results in an estimate (still assuming an isother- 
mal sphere) of the minimum mass of a star that can reach the cluster core 
in time t due to dynamical friction: 

Thus, at time t and for a given mass to, there is a maximum radius R 
inside of which stars of mass to will have segregated to the core. The 
stars contributing to the growth of the runaway are likely to be among 
those more massive than m&f, because their number density in the core is 
enhanced by mass segregation, their collision cross sections are larger, and 
they contribute more to <5m co n when they do collide. 

The shape of the central mass function of a segregated cluster is not 
trivial to derive 12 . In thermal equilibrium, the central number densities of 
stars of different masses would be expected to scale as 

n (m) ~ to 3/2 — , (1.34) 

dm 

where dN /dm is the global initial mass function, which scales roughly as 
m~ 27 at the high-mass end (to £ 10M Q ). The distribution of secondary 
masses (i.e. the masses of the lighter stars participating in collisions) does 
not follow the above simple relation. Rather, we find that stars in the 
core do not reach thermal equilibrium (a result generally consistent with 
earlier findings by Chernoff and Weinberg 1990 and Joshi, Nave & Rasio 

11 The author realizes that the simple mention of primordial binaries opens-up an entire 
discussion which would require several pages, which I try to prevent. 

12 In that case, the mass function in the core takes on a rather curious form: the mass 
functions for stars with masses m ^ m^t and m ^ m^f have roughly the same slopes 
as the initial mass function, but the more massive stars are overabundant because they 
have accumulated in the cluster center. 
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2001), and that the dynamical nature of the collisional processes involved 
means that more massive stars tend to be consumed before lower-mass stars 
arrive in the core. In addition, most collisions involve three-body binaries 
and interactions with higher order multiples in a multi-mass environment. 
Empirically, we find that the secondary mass distribution is quite well 
fit by a power-law, dN/dm ex m~ 2 ' 3 (coincidentally very close to a Salpeter 
distribution). Integrating this expression from a minimum mass of maf 
(and ignoring the upper limit) results in a mean mass increase per collision 
of 

£ro co u — 4m df . (1.35) 

If we neglect stars with masses less than m<if and substitute Eg . 11. 51 into 
Eg . 11. 331 and Eo. ll.35l then the mass increase per collision can be written as 

(5m coll ~4— (m)lnA. (1.36) 

This quantity remains rather constant over the entire collisional lifetime, 
e.g., about 3Myr (see also § ll.3.3.2|) . 

1.3.1.3 Lifetime of a cluster in a static tidal field 

With simple expressions for 7V co ii and 8m co \\ now in hand, we return to 
the determination of the runaway growth rate fEo. ll.3()j) . The evaporation 
of a star cluster which fills its Jacobi surface in an external potential is 
driven by tidal stripping. Portegies Zwart et al (2001a) have studied the 
evolution of young compact star clusters within ~ 200 pc of the Galactic 
center. Their calculations employed direct A-body integration, including 
the effects of both stellar and binary evolution and the (static) external 
influence of the Galaxy, and made extensive use of the GRAPE-4. They 
found that the mass of a typical model cluster decreased almost linearly 
with time: 

M = M (l - -^-] . (1.37) 

\ tdisr / 

Here Mq is the mass of the cluster at birth and idisr is the cluster's disrup- 
tion time. Portegies Zwart et al. (2001a) found that their model clusters 
dissolved within about 30% of the two-body relaxation time at the tidal 
radius (defined by substituting the tidal radius instead of the virial radius 
in Ea . 11. 5(1 . In terms of the half- mass relaxation time, this translates to 
£disr = 1.6-5.4< r t, depending on the initial density profile (the range cor- 
responds to King [1966] dimensionless depths Wo = 3-7; more centrally 
condensed clusters live longer). 

Substituting Eos. 02l and 06l into Eo. OUl and defining M = N (m) 
to rewrite Eo. 11.371 in terms of the number of stars in the cluster, we find 

dM run 3 JV(m)lnA 

dt Jc t 
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= 4 x 10- 3 /cM lnA (~ - -?-) . (1.38) 

V * 'disr / 



Integrating from t — t cc to t — t^isr results in 

M run = m sccd + 4 x 10- 3 / c A/ ln A [in [^] + -^- 



(1.39) 



Here m see d is the seed mass of the star which initiates the runaway growth, 
most likely one of the most massive stars initially in the cluster. With 
i cc ~ 0.2i rt , E all. 391 reduces to 

M run = m socd + 4 x 10" 3 /cMo«:ln A , (1.40) 

where k ~ ln (£ disr /£ cc ) + t cc /t disr - 1 ~ 1. 

In figur cll~T7l we present this relation in the form of a solid curve. We 
comment further on the left side of this figure where we extrapolate the 
relation in Ea. ll.4Ul to galactic nuclei masses. 

The maximum mass of the runaway merger for clusters which are 
disrupted by inspiral (which of course always destroys the cluster before it 
reaches the center) may be calculated by replacing £ d ; sr in Ea . 11.391 bv T{. 
The right-hand side of that equation then becomes a function of 

t„ U<W \ R ) V M ) ■ ' 

We can also estimate the maximum initial distance from the Galactic 
center for which core collapse occurs (and hence runaway merging may 
begin) before the cluster disrupts by setting T{ — t cc . The result is 
Ri > 0.0025pc(i?A//[pcM Q ]) a71 . For R = 0.25 pc and M = 10 5 M Q , we 
find Ri £ 3.3 pc. 

1.3.2 Calibration with AT-body simulations 

The development of the GRAPE (see S ll.1.4.41) family of special-purpose 
computers makes it relatively straightforward to test and tune the above 
theory using direct A-body calculations. The N-hody calculations per- 
formed by Portegies Zwart & McMillan (2002) span a broad range of initial 
conditions in the relevant part of parameter space. The number of stars 
varied from Ik (1024) to 64k (65536). The initial conditions explored by 
Portegies Zwart et al (2004) ranges from 131,072 to 585,000 stars. Giirkan 
et al (2004) performed similar calculations using a Monte-Carlo A-body 
code, they adopted a considerably larger number of particles. 

Initial density profiles and velocity dispersion for the models were 
taken from Heggie-Ramamani models (Heggie & Ramamani 1995) with 
Wq ranging from 1 to 7, and from King (1966) models with Wo = 3 
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Figure 1.17. The mass after a period of runaway growth as a function of the 
mass of the star cluster. The solid line is M Tun = 30 + 8 x 10 _4 MolnA (see 
Eg . 11. 401 with f c — 0.2, 7 = 1 and lnA= InMo/M©, where M is the initial mass 
of the cluster or 1O 6 M0, whatever is smaller). This relation may remain valid 
for larger systems built up from many clusters having masses ,$ 10 Mq. For 
clusters with Mo ^ 10 M© we therefore extend the relation as a dashed line. The 
logarithmic factor, however, remains constant, as it refers to the clusters out of 
which the bulge formed, not the bulge itself (Ebisuzaki et al 2001). The bottom 
dashed line shows 0.01M run . The five error bars to the left give a summary of 
the results presented by Portegies Zwart & McMillan (2002), the two right most 
error bars are taken from Portegies Zwart et al (2004). The error bar indicating 
MGG11 gives the result of estimates for the mass of the intermediate mass black 
hole in the M82 star cluster MGG11 (Matsumoto & Tsuru 1999; McCrady et al 
2003). The Milky Way is represented by the asterisk using the bulge mass from 
Dwek (et al. 1995) and the black hole mass from Eckart & Genzel (1997) and 
Ghez (2000). Bullets and triangles (upper right) represent the bulge masses and 
measured black hole mass of Seyfert galaxies and Quasars, respectively (both 
from Wandel 1999; 2001). The dotted lines gives the range in solutions to a least 
squares fits to the bullets and triangles (Wandel 2001). 



15. At birth, the Heggie-Ramamani clusters were assumed to fill their 
zero- velocity ( Jacobi) surfaces in the Galactic tidal field, while the classical 
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King models were assumed to be isolated. In most cases we adopted an 
initial mass function between 0.1 M Q and 100 M suggested for the So- 
lar neighborhood by Scalo (1986) and Kroupa, Tout & Gilmore (1990). 
However, several calculations were performed using power-law initial mass 
functions with exponents of -2 or -2.35 (Salpeter) and lower mass limits of 
1M . 

1.3.2.1 Collision rate during phase A 

In all calculations, the first collision occurred shortly after the formation 
of the first £ 10 kT binary by a three-body encounter, i.e. close to the 
time of core collapse. When stars were given unrealistically large radii (100 
times larger than normal), the first collisions occurred only slightly (about 
5%) earlier. 

As discussed earlier, the first star to experience a collision was generally 
one of the most massive stars in the cluster; this star then became the target 
for further collisions. In models where the core collapse time exceeds about 
3 Myr the target star explodes in a supernova before experiencing runaway 
growth. The collision rates in these clusters were considerably smaller 
than for clusters with smaller relaxation times (see Fig. ll,18|) . As discussed 
in more detail in §4, the onset of stellar evolution terminates the collision 
process; premature disruption of the cluster also ends the period of runaway 
growth. 

The first physical collisions occur at the moment of core collapse. The 
cluster then enters a phase which is dominated by stellar collisions. In 
particular one single object experiences many repeated collisions, giving rise 
to a collision runaway. We identify this particular object as the designated 
target. In our models collisions tend to increase the mass of the collision 
product; only a small fraction of the mass of the incoming star is lost; the 
designated target therefore tends to increase in mass. 

The number of collisions in the direct ./V-body simulations ranged from 
to 100. Fig. ll.18l shows the mean collision rate 7V co ii per star per million 
years as a function of the initial half-mass relaxation time. The solid line 
in Fig. 11.181 is a fit to the simulation data, and has 

Kou = 2.2 x 10- 4 — , (1.42) 

trt 

for i rt ^S 20 — 30 Myr, consistent with our earlier estimate fEq. ll.32JI if 
/ c = 0.2. The quality of the fit in Fig. ll.18l is quite striking, especially 
when one bears in mind the rather large spread in initial conditions for 
the various models. See however the prominent square to the right, which 
is about a factor ~ 3 above the fitted curve. This discrepancy is mainly 
caused by the high average stellar mass in these models and by the use of 
much more concentrated King models. 
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Figure 1.18. Mean collision rate f co ii = Ncoil/Ntiast as function of initial relax- 
ation time for the simulations performed by Portegies Zwart et al. 1999; 2004 
and Portegies Zwart & McMillan 2002). Here £i as t is the time of the last collision 
in the cluster. The open circles give the results of systems which are isolated 
from the Galactic potential (see Portegies Zwart et al 1999). Vertical bars rep- 
resent Poissonian 1-cr errors. The solid line is a least squares fit to the data (see 
Eq. ll.42t . The strong reduction in the collision rate for cluster with an initial 
relaxation time t It ^ 30Myr is probably real. Note however the filled square 
to the right from the calculations of Portegies Zwart et al (2004), which gives a 
higher collision rate due to the high initial concentration of the models and the 
top-heavy mass function. 
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The collision runaway phase lasts until about 3.3 Myr, at which time 
the designated target tends to collapse to a black hole. 

In fig. II. 191 we present the cumulative distribution of the number of 
mergers of some of the calculations by Portegies Zwart et al (2004). 
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Figure 1.19. Cumulative distributions of the number of mergers which occurred 
within time t. Panel a) gives the results of the calculations performed with 
starlab, panel b) is from NB0DY4. The dotted, solid and dashed curves are for 
the models with Wo = 8 (in the right panel we show the results for Wo — 7), 
Wo = 9 and Wo = 12. The dash-3-dotted curve in panel a) is for a model with 
Wo — 12 with 10% primordial binaries. The dash-3-dotted curve in panel b) is 
for a model with Wo = 9 using 585 000 stars, distributed according to a Kroupa 
(2001) initial mass function between O.lMoand 100 Mq. 

Figur eTT.201 shows the cumulative mass distributions of the primary 
(more massive) and secondary (less massive) stars participating in colli- 
sions. We include only events in which the secondary experienced its first 
collision (that is, we omit secondaries which were themselves collision prod- 
ucts). In addition, we distinguish between collisions early in the evolution 
of the cluster and those that happened later by subdividing our data based 
on the ratio r = £ C oiiAdf, where t co \\ is the time at which a collision oc- 
curred and idf is the dynamical friction time scale of the secondary star 
(see Ea. ll.21|l . The solid lines in Figur eTOni show cuts in the secondary 
masses at r ;$ 1, r ^E 5 and r < oo (rightmost line). The mean secondary 
masses are (to) = 4.0 ± 4.8 M Q , 8.2 ± 6.5 and (to) = 13.5 ± 8.8 M Q for 
t ,$ 1, 5 and oo, respectively. 

The distribution of primary masses in Figure lT.2UI (dashed line) hardly 
changes as we vary the selection on r. We therefore show only the full 
(r ;$ oo) data set for the primaries. In contrast, the distribution of sec- 
ondary masses changes considerably with increasing r. For small r, sec- 
ondaries are drawn primarily from low-mass stars. As r increases, the 
secondary distribution shifts to higher masses while the low-mass part 
of the distribution remains largely unchanged. The shift from low-mass 
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( ;$ 8M ) to high-mass collision secondaries ( £ 8M Q ) occurs between 
r = 1 and t = 5. This is consistent with the theoretical arguments pre- 
sented in Sec. ll.3.1~21 During the early evolution of the cluster (r ;$ 1), 
collision partners are selected more or less randomly from the available 
(initial) population in the cluster core; at later times, most secondaries are 
drawn from the mass-segregated population. 

Interestingly, although hard to see in Fig. 11.201 all the curves are well 
fit by power laws between ~ 8M and ~ 80M Q (0.8M o and 3OM for 
the leftmost curve). The power-law exponents are —0.4, —0.5 and —2.3 for 
t £ 1, 5, and oo, and —0.3 for the primary (dashed) curve. (Note that the 
Salpeter mass function has exponent —2.35.) 
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Figure 1.20. Cumulative mass distributions of primary (dashed line) and sec- 
ondary (solid lines) stars involved in collisions. Only those secondaries experienc- 
ing their first collision are included. From left to right, the solid lines represent 
secondary stars for which r = t co \i/t&t 
included in each curve are 18 (for t co n/tdf 
The dotted line gives a power-law fit with the Salpeter exponent (between 5 Mq 
and 100 Mq) to the right most solid curve (r ,$ oo). 



$, 1, 5, and oo. The numbers of collisions 
$ 1), 42 and 95 (rightmost two curves). 



Figur elT. 1 7l shows the maximum mass reached by the runaway collision 
product as function of the initial mass of the star cluster. Only the left side 
(log M/M ;$ 7) of the figure is relevant here; we discuss the extrapolation 
to larger masses in Sec. ll.4.2l The N-body results are consistent with the 
theoretical model presented in Eg . 11. 391 
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1.3.3 Simulating the star cluster MGG11 

In a recent publication (Portegies Zwart et al, 2004) we simulate a well 
observed star cluster in the starburst Galaxy M82. Interestingly enough, 
this cluster has experienced a prominent phase A evolution, and is currently 
in a phase B. In this § I will report some of the interesting results about 
these simulations. Note that we have been quoting these results in earlier 
occasions in the text, but here we review the global results. 

The observed mass function, as reported by McCrady et al (2003) for 
the M82 star cluster MGG11 is consistent with a Salpeter power-law (with 
x = —2.35) between 1M and an upper limit which corresponds to the 
age of the cluster of 7 to 12Myr. By that time all stars more massive 
than 17-25 M Q have experienced a supernova. For the IMF we adopt the 
same Salpeter slope and lower mass limit of 1 M , but extend the upper 
limit to 100 M Q . This IMF has an average mass of (to) ~ 3.1 M . If we 
assume that at an age of 12-7 Myr all stars between 17-25 M and 100 M Q 
are lost from the cluster by supernova explosions, the mean stellar mass 
drops to (to) ~ 2.4 - 2.6 M . With a current total mass of 3.5 x 10 5 M 
the cluster would contain 130 000 to 140 000 stars. For clarity we decided 
to select 128k (131072) single stars, resulting in an initial mass of about 
406,000 M . 

McCrady et al (2003) measured the projected half light radius of the 
cluster, rhp = 1.2 pc. De-projection of the half light radius depends on the 
density profile. For King (1966) models in the range Wo — 3 — 12 it turns 
out that rjj = 0.75 — 0.88rh p (which is somewhat larger than Spitzer's, 1987 
r h ~ f r hp)- A number of initial test simulations indicate that over a time 
span of 7-12 Myr the projected half light radius of the selected IMF and 
number of stars, the cluster expands by about a factor 1.3. We then adopt 
an initial half mass radius for our models cluster of r^m = 1.2 pc. 

We ignore the effect of the tidal field of M82. The star cluster is 
located at a distance of about 160 pc from the dynamical center of the 
Galaxy, assuming a distance of 3.6 Mpc (Freedman et al. 1994) 13 . With 
the relatively small mass of M82 of ~ 1O 9 M the gravitational force of 
the Galaxy is negligible compared to the self gravity of the stars within the 
cluster. 

We performed several calculations starting from King models with 
different central concentrations Wo- We also performed one run with 10 
per cent primordial binaries and one run starting with 585 000 stars and a 
Salpeter IMF between 0.1 M and 100 M . 

To qualify the results we make the distinction between clusters with 
a high central concentration and clusters with low concentration. Since 
the initial half mass radius is the same for all models, we varied the den- 

13 Freedman et al. measured a distance of 3.63 ± 0.34 Mpc to M81, the neighboring 
galaxy of M82, using Cepheids. 
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Table 1.4. Resulting collisions for our model clusters as a function of the initial 
density profile expressed in the King parameter Wo (identified on the first column) 
and the average core density is presented in the second column. The third column 
gives the number of collisions which occurred in 12 Myr. The fourth and fifth 
columns give the number of collisions in which one particular star participates 
and the average mass of the star with which it collides. The last three columns 
give the number of collision and the mean mass of the more massive and lower 
mass stars participating in these collisions. 
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sity profile. For the density profile we adopted King (1966) models. We 
draw the empirical distinction between high concentration models having 
Wo ^ 9; low concentration models have Wq ^S 8. With the adopted half- 
mass radius and total mass the high concentration cluster Wo ^ 9 models 
have core density logp c > 6M pc~ 3 . 

In figur eTl"2"Tl we present the evolution of the central density of various 
starlab models with Wq — 7, 8 and 9. As discussed, in the low con- 
centration models core collapse is arrested by the copious mass loss from 
the evolving massive stars. In the high concentration clusters core collapse 
occurs early enough that the process is little affected by stellar mass loss. 



1.3.3.1 Clusters with logp c > 6M pc~ 3 

According to Portegies Zwart et al (1999), who studied similar clusters as 
Portegies Zwart & McMillan (2002), the high collision rate is mainly the 
result of binaries created in 3-body encounters during core collapse. In 
our latest models, however, only about half (0.52 ± 0.1) the collisions are 
preceded by the formation of a binary, the other half are direct hits in which 
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Figure 1.21. Evolution of the central density for star clusters with Wo — 7 
(dash-3-dotted line), W = 8 (dahses), W = 9 (dotted curve). The W = 9 
model shows a clear core collapse at about 0.7 Myr, indicated with the arrow. 
The less concentrated models (Wo — 8 and Wo — 7) show a very shallow start 
of a core collapse near t ~ 3 Myr, but in these cases the collapse of the core is 
arrested by stellar mass loss. In fig. ll.13l we present the evolution of the core 
radius for the same calculations. 



no third star was bound to the two stars which ultimately collided. Though 
seemingly a detail, it has far reaching consequences for the interpretation 
of the collisional growth which plays an important role in the evolution of 
these systems. 

After the first epoch of runaway growth and the collapse of the des- 
ignated target, an epoch of about 3 Myr starts in which the collision rate 
drops dramatically. This phase is visible in fig. ll.lijl between ~ 3.3 Myr and 
- 6 Myr. 

This quiet phase lasts until the first M5I hyper giants appear, which, 
in our stellar evolution model, happens at a turn-off mass of ~ 25 M 
(at ~ 6 Myr) . (Note that in our stellar evolution model SeBa stars below 
~ 25 M Q turn into neutron stars where more massive stars become black 
holes, see Portegies Zwart & Verbunt 1996) By this time the collision rate 
picks up again to 2-3 per Myr. The spectral type MOI - M5I star dominate 
the collision rate; the designated target, now an intermediate mass black 
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hole, participates in only about one-third of the collisions. This phase lasts 
until about 9-12 Myr, after which the rate drops below 0.3 collisions per 

Myr. 

1.3.3.2 Clusters with logp c < 6M pc~ 3 

In low density cluster Wq £ 8 the initial phase of runaway growth is absent. 
The subsequent collisional phase at an age £ 6 Myr, however, occurs at 
a comparable rate as in the high density simulations (see fig. ll.20|l . As 
discussed in the previous § there is no designated target in this phase; all 
collisions tend to happen between massive stars (or stellar remnants) of 
which at least one component is evolved (spectral type MOI to M5I). At 
a rate which becomes gradually smaller for less concentrated clusters: 3.4 
collisions per Myr for W = 9, 3.2 for W = 8, 1.6 for W = 7, 1.5 for 
Wo = 5 and 0.5 for Wo = 3. Curiously enough the Wo = 12 models have a 
collision rate of only 2.3 per Myr, which is lower than the less concentrated 
W Q — 9 models. 

The average star with which a collision occurs (counting the least 
massive of the two stars) depends on the initial cluster density profile, 
ranging from a 5.3 M Q for W = 5, 6.2 M Q for W = 7, 9.7 M Q for W = 8 
and ~ 30 for both Wq — 9 and Wq = 12 (consistent with the expectations 
of Portcgies Zwart & McMillan 2002) 

1.3.4 Black hole ejection in phase B and C cluster with i rt ^ 100 Myr 

Here we discuss the evolution of star clusters in which the early phase 
is dominated by stellar mass loss and the subsequent evolution by the 
stellar interaction; tse < i~sd- In this regime or parameter space, the most 
massive stars evolve before the structure of the star cluster has appreciably 
changed, i.e; no intermediate mass black hole can form via the scenario 
discussed in § 11.3.11 The consequence is that the most massive stars turn 
into black holes and neutron stars before they had a chance to sink to the 
cluster center by dynamical friction. This regime is valid for most globular 
clusters, and possibly many open star clusters. 

Upon the birth of a cluster we assume that the stars populate the 
initial mass function from the hydrogen burning limit (~ 0.1 Mg) all the 
way to the most massive stars currently observed in the Galaxy, which is 
about 100 Mg. Stars with Solar abundance between 50 M and 100 M Q 
leave the main sequence at an age of about 3.7 Myr to 3.3 Myr, to explode 
in a supernova a few hundred thousand years later. The total mass in this 
range is about 1%, and the cluster therefore loses approximately 1% of its 
total mass in less than half a million years. 14 

14 For an entire star cluster a 1% mass loss is not very dramatic, and simply causes the 
cluster to expand by about the same fraction. For a cluster core, which contains less 
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The first black holes are produced at about the same time. Black hole 
formation proceeds to about 7-9 Myr, until all stars with initial masses 
exceeding 20-25 M (Maeder 1992; Portegies Zwart et al. 1997). have 
collapsed to black holes. Assuming a Scalo (1986) mass distribution with a 
lower mass limit of 0.1 M Q and an upper limit of 100 M© about 0.071% of 
the stars are more massive than 20 Mg, and 0.045% are more massive than 
25 Mq. A star cluster containing N stars thus produces ~6x 10 _4 AT black 
holes. Known Galactic black holes have masses TObn between 6M© and 
18 Mq (Timmes et al. 1996). For definiteness, we adopt rribh = 10 Mq. 

A black hole is formed in a supernova explosion. If the progenitor is 
a single star (i.e. not a member of a binary), the black hole experiences 
little or no recoil and remains a member of the parent cluster (White & 
van Paradijs 1996). 15 If the progenitor is a member of a binary, mass loss 
during the supernova may eject the binary from the cluster potential via 
the Blaauw mechanism (Blaauw 1961), where conservation of momentum 
causes recoil in a binary as it loses mass impulsively from one component. 
The Blaauw velocity kick can be as large as the relative orbital velocity 
of the pre-supernova binary (see e.g. Nelemans et al (2001). The escape 
speed is ~ 40 km s _ for a young globular cluster, and somewhat smaller 
for YoDeCs (see Tab ll.lll . Such high recoil velocities are generally achieved 
only if the binary loses ^S 50% of its total mass in the supernova, and 
if its orbital period is initially quite small ( ,$ 2yr) (Portegies Zwart & 
Verbunt 1996). The binary frequency in globular clusters is between 5 
and 40% (Rubenstein 1997), and less than 30% of binaries have orbital 
periods smaller than 2 years (Rubenstein & Bailyn 1997; 1999); we assume 
the same distributions of orbital parameters for binaries in YoDeCs. We 
therefore estimate that no more than ~ 10% of black holes are ejected 
from the cluster immediately following their formation (i.e. the black hole 
retention fraction is Z, 90%). Note the remarkable difference here with the 
retention fraction of neutron stars, which is less than about 10% (Drukier 
1996). 

After ~ 40 Myr the last supernova has occurred, the mean mass of 
the cluster stars is (m) ~ 0.56 M (Scalo 1986), and black holes are by 
far the most massive objects in the system. Mass segregation cause the 
black holes to sink to the cluster core in a fraction ~ (to) /rribh of the half- 
mass relaxation time scale (see Ea. ll.21jl . For a young populous cluster, 
the relaxation time is on the order of 10 Myr; for a globular cluster it is 
about lGyr (see Tab. Ojl , 

By the time of the last supernova, stellar mass loss has also significantly 
diminished and the cluster core starts to contract, enhancing the formation 

than 5% of the total cluster mass (for Wo ^i 8), a 1% mass loss may drive a substantial 
expansion of the core. 

15 There are some arguments which indicate that black hole may received small 'impulse' 
kicks, relative to neutron stars. 
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of binaries by three-body interactions. Single black holes form binaries 
preferentially with other black holes (Kulkarni et al. 1993), while black 
holes born in binaries with a lower-mass stellar companion rapidly exchange 
the companion for another black hole. The result in all cases is a growing 
black-hole binary population in the cluster core. Once formed, the black- 
hole binaries become more tightly bound through superelastic encounters 
with other cluster members (Heggie 1975; Kulkarni et al. 1993; Sigurdsson 
& Hernquist 1993). On average, following each close binary-single black 
hole encounter, the binding energy of the binary increases by about 20% 
(Hut et al. 1992); roughly one third of this energy goes into binary recoil, 
assuming equal mass stars. The minimum binding energy of an escaping 
black-hole binary may then be estimated as 

£6, mln ~ 36 Wo^feT, (1.43) 

(to) 

where |fcT is the mean stellar kinetic energy and Wq — (m)\(f)o\/kT is 
the dimensionless central potential of the cluster (King 1966). By the 
time the black holes are ejected, (to) ~ 0.4 M Q . Taking Wq ~ 10 we find 
£&,min~ 10000 kT. 

Portegies Zwart & McMillan (2000) have tested and refined the above 
estimates by performing a series of N-hody simulations within the Starlab 
software environment using the special-purpose computer GRAPE-4 to 
speed up the calculations (see S I1.1.4.4J I. For most (seven) of these calcula- 
tions we used 2048 equal-mass stars with 1% of them ten times more mas- 
sive than the average; two calculations were performed with 4096 stars. One 
of the 4096-particle runs contained 0.5% black holes; the smaller black-hole 
fraction did not result in significantly different behavior. They also tested 
alternative initial configurations, starting some models with the black holes 
in primordial binaries with other black holes, or in primordial binaries with 
lower-mass stars. 

The results of these simulations may be summarized as follows. Of 
a total of 204 black holes, 62 (~ 30%) were ejected from the model clus- 
ters in the form of black-hole binaries. A total of 124 (~ 61%) black 
holes were ejected single, and one escaping black hole had a low-mass 
star as a companion. The remaining 17 (~ 8%) black holes were re- 
tained by their parent clusters. The binding energies E\, of the ejected 
black-hole binaries ranged from about 1000 kT to 10000 kT in a distribu- 
tion more or less flat in log Ei,, consistent with the assumptions made by 
Hut et al. (1992). The eccentricities e followed a roughly thermal distri- 
bution [p(e) ~ 2e], with high eccentricities slightly overrepresented; the 
mean eccentricity was (e) = 0.69 ± 0.10. The 17 binaries with the low- 
est binding energies (log 10 Ej, < 3.5) had on average higher eccentrici- 
ties ((e) = 0.78 ± 0.05) than the more tightly bound binaries ((e) = 
0.62 ± 0.11). About half of the black holes were ejected while the par- 
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ent cluster still retained more than 90% of its birth mass, and /t 90% of 
the black holes were ejected before the cluster had lost 30% of its initial 
mass. These findings are in good agreement with previous estimates that 
black-hole binaries are ejected within a few Gyr, well before core collapse 
occurs (Kulkarni et al. 1993; Sigurdsson & Hernquist 1993). Recently 
Merritt et al. (2004) confirmed these findings with their own simula- 
tions using NB0DY6++, (Hcmsendorf et al 2002; and can be downloaded 
from ftp : //ftp . ari . uni-heidelberg . de/pub/staf f /spurzem/nb6mpi/ 
a parallel version of NB0DY6 (Aarseth 1999). 

Portegies Zwart & McMillan (2000) have performed additional cal- 
culations incorporating a realistic (Scalo 1986) mass function, the effects 
of stellar evolution, and the gravitational influence of the Galaxy. These 
model clusters generally dissolved rather quickly (within a few hundred 
Myr) in the Galactic tidal field. We found that clusters which dissolved 
within ~ 40 Myr (before the last supernova) had no time to eject their black 
holes; however, those that survived for longer than this time were generally 
able to eject at least one close black-hole binary before dissolution. 

Based on these considerations, we conservatively estimate the number 
of ejected black-hole binaries to be about 10 -4 iV per star cluster, more or 
less independent of the cluster lifetime. 

1.3.4-1 Characteristics of the black-hole binary population 

The energy of an ejected binary and its orbital separation are coupled to 
the dynamical characteristics of the star cluster. For a cluster in virial 
equilibrium, we have 
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In computing the properties of the black-hole binaries resulting from 
cluster evolution, it is convenient to distinguish three broad categories of 
dense stellar systems: (1) YoDeCs, (2) globular clusters, and (3) galactic 
nuclei. Tabl elLTI lists characteristic parameters for each category. The 
masses and virial radii of globular clusters are assumed to be distributed 
as independent Gaussian with means and dispersions as presented in the 
table; this assumption is supported by correlation studies (Djorgovski & 
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Meylan 1994). Table H~Tl also gives estimates of the parameters of globular 
clusters at birth (bottom row), based on a recent parameter-space survey 
of cluster initial conditions (Takahashi & Portegies Zwart 2000; Baumgardt 
& Makino 2003); globular clusters which have survived for a Hubble time 
have lost £ 60% of their initial mass and have expanded by about a factor 
of three. We draw no distinction between core-collapsed globular clusters 
(about 20% of the current population) and non-collapsed globulars — the 
present dynamical state of a cluster has little bearing on how black-hole 
binaries were formed and ejected during the first few Gyr of the cluster's 
life. 

The above described process causes globular clusters to be depleted of 
black holes. At most one binary containing two black holes can be present 
in any globular cluster. In S ll.4.4l we will further discuss the consequences 
of the ejected black holes and black hole binaries, as the latter may become 
important sources for gravitational wave detectors. 16 

1.4 Discussion and further speculations 

Now we have discussed the basic principles of star cluster dynamics with 
black holes. In this section we further discuss some of the consequences 
and observables to which this theory can be tested. 

1.4.1 Turning an intermediate mass black hole in an X-ray source 

An IMBH in a stellar cluster with velocity dispersion a dominates the po- 
tential well within its radius of influence i?; = Gm^/a 2 ; inside i?i the 
orbits are approximately Keplerian, and the stars are distributed accord- 
ing to a power law (see Ea. ll.2.l"Tl and also (Bahcall & Wolf) . N-body 
calculations show that a = 1.5 (Baumgardt et al 2004), and we assume 
this value in the following. 

Stars can reach an orbit with peribothron of order of the tidal radius 
by energy diffusion or by angular momentum diffusion. However, stellar 
collisions become the dominant dynamical process within the collision ra- 
dius r co ii ~ (Wbh/w)r (Frank & Rees 1976), disrupting stars within that 
region, and making energy diffusion within r co ii implausible. For this rea- 
son we concentrate on tidal capture of stars on very eccentric orbits. When 
the star arrives at peribothron a certain amount of energy AEt is invested 
in raising tides, causing the orbit to circularize. It is hard to know how 
the star dissipates the tides after the repetitive encounters with the IMBH. 
Two extreme models of "squeezars" (stars that are "squeezed" by the tidal 
field) are studied by Alexander & Morris (2003), namely "cold squeezars", 

16 Today's globular clusters may undergo a similar ejection phase, which can be seen in 
the high proportion of recycled pulsars and the interestingly possibility of a high merger 
rate of white dwarfs (Shara & Hurley, 2002). 
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which puff up, or "hot squeezars" , which are heated only in their outer 
layers and radiate their excess energy effectively. 

Hopman et al (2003) calculated that the probability that an intermedi- 
ate mass black hole captures a stellar companion via this process (assuming 
the "hot squeezar" model) has a reasonable probability, making it likely 
that we could in fact observe one or two of such objects. Once captured 
the companion star is likely to fill its Roche-lobe while still on the main- 
sequence. 

From the moment of first Roche-lobe contact the evolution of the bi- 
nary is further determined by mass transfer from a lower mass (secondary) 
star to the much more massive (primary) IMBH. This process is driven 
by the thermal expansion of the donor and the loss of angular momentum 
from the binary system. Mass transfer then implies that the donor fills its 
Roche lobe (i?don = Rri) and continues to do so, i.e., i?don = Rm- We 
also assume that, as long as Ldisc < -^Edd all the mass lost from the donor 
via the Roche lobe is accreted by the black hole, M = —to. If the mass 
transfer rate exceeds the Eddington limit, the remaining mass is lost from 
the binary system with the angular momentum of the accreting star. 

We estimate the X-ray Luminosity during mass transfer using the sim- 
ple model of Kording, Falcke & Markoff (2002). They argue that in the 
hard state (m > rh C rit) L x ~ O.lmc 2 , where at lower accretion rates the 
source becomes transient (see also Kalogera et al, 2004) with an average 
luminosity of L x ~ 0.1c 2 m /m C rit- For m cr ; t we use the relation derived 
by Dubus et al (1999) 

1.4.2 Speculation on the formation of supermassive black holes 

A million-solar- mass star cluster formed at a distance of ;$ 30 pc from 
the Galactic center can spiral into the Galactic center by dynamical fric- 
tion before being disrupted by the tidal field of the Galaxy (see Gerhard 
2001). Only the densest star clusters survive to reach the center. These 
clusters are prone to runaway growth and produce massive compact objects 
at their centers. Upon arrival at the Galactic center, the star cluster dis- 
solves, depositing its central black hole there. Black holes from in-spiraling 
star clusters may subsequently merge to form a supermassive black hole. 
Ebisuzaki et al. (2001) have proposed that such a scenario might explain 
the presence of the central black hole in the Milky way galaxy. 

If we simply assume that bulges and central supermassive black holes 
are formed from disrupted star clusters, this model predicts a relation be- 
tween black hole and bulge masses in galaxies similar to the expression 
(Eq. ll.4()|l connecting the mass of an IMBH to that of its parent cluster. 
However, the ratio of stellar mass to black-hole mass might be expected 
to be smaller for galactic bulges than for star clusters, because not all star 
clusters produce a black hole and not all star clusters survive until the 



Discussion and further speculations 




53 

iM/dt 

-5 



0.5 1 

log(t/Myr) 



Figure 1.22. Estimated X-ray luminosity (log erg/s) as a function of time 
for a 1000 Mq accreting black hole. The solid, dashed and dotted curves are 
for a 15 Mq , 10 M© and a 5 M© donor which started Roche-lobe overflow at the 
zero-age main-sequence. (Figure from Hopman, Portegies Zwart & Alexander, 
2004) 



maximum black hole mass is reached. We would expect, however, that the 
general relation between the black hole mass and that of the bulge remains 
valid. 

Figure H~T7l shows the relation between the black hole mass and the 
bulge mass for several Seyfert galaxies and quasars. The expression derived 
in Sec. II. 3. PI and the results of these TV-body calculations (Sec. 1 1.3. "2]) are 
also indicated. The solid and dashed lines fEa. ll.40]l fit the A^-body calcula- 
tions and enclose the area of the measured black hole mass-bulge masses. 
On the way, the solid curve passes though the two black-hole mass esti- 
mates, for M82 star cluster MGG-11. We note that the observed relation 
between bulge and black hole masses has a spread of two orders of magni- 
tude. If this bold extrapolation really docs reflect the formation process of 
bulges and central black holes, this spread could be interpreted as a vari- 
ation in the efficiency of the runaway merger process. In that case, only 
about one in a hundred star clusters reaches the galactic center, where its 
black hole is deposited. 

The dashed curve is an extrapolation beyond 10 7 M the range where 
we think that Ea. ll.40l is applicable. There appears to be a very interesting 
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relation for galactic nuclei between the velocity dispersion and the mass 
of the central black hole (Ferrarese & Merritt 2000; Gebhardt et al 2000, 
2001); and Merritt & Ferrarese (2001). I do not wish to claim that the 
process described in Ea. ll.40l and the formation of supermassive black holes 
in galactic nuclei (see however Ebisuzaki et al 2001), but just point out how 
well the same theory explains the relation between nuclear black hole mass 
and the velocity dispersion in the nucleus. At the moment of writing there 
more than a dozen alternative explanations for this interesting relation 
(see for example Haehnelt et al. 2000; Dokushaev et al 2002; Boroson et 
al 2002); Green et al (2004);), and this number is growing almost daily 17 . 

1.4.3 Is the globular cluster M15 a special case? 

The presence of an intermediate mass black hole in the core collapsed glob- 
ular cluster M15 has been vividly debated over the last several decades 
(Illingworth & King, 1977a, 1977b) and is still discussed to the present 
day (van der Marel et al 2002). If M15 contains an intermediate mass 
black hole, as has recently been suggested by some observers (Gerrisen et 
al, 2002; later retracted and subsequently corrected in Gerrisen 2003) and 
argued against by theorists (Baumgardt et al 2003), it is unlikely to have 
formed via the scenario discussed in § 11.3.11 The cluster's initial relaxation 
time probably exceeded our upper limit of ~ 100 Myr. The current half- 
mass relaxation time of M15 is about 2.5 Gyr (Harris 1996), which is far 
more than our 100 Myr limit for forming a massive central object from a 
collision runaway (see 8 11.3.11 for details). 

An intermediate mass black hole in the globular cluster M15 may have 
been formed by a different scenario. An alternative is provided by Miller & 
Hamilton (2001), who describe the formation of massive (~ 10 3 M ) black 
holes in star clusters with relatively long relaxation times. In their model 
the black hole grows very slowly over a Hubble time via occasional collisions 
with other black holes, in contrast to the model described here, in which 
the runaway grows much more rapidly, reaching a characteristic mass of 
about 0.1% of the total birth mass of the cluster within a few megayears. 

One possible way around M15's long relaxation time may involve the 
cluster's rotation. Gebhardt (2000; 2001; private communication) has mea- 
sured radial velocities of individual stars in the crowded central field, down 
to two arcsec of the cluster center. He finds that, both in the central part 
of the cluster (R < 0.1i?hm) and outside the half mass radius, the average 
rotation velocity is substantial (i> ro t ~ 0.5(v)). Rotation is quickly lost in 
a cluster (Baumgardt et al 2003), so to explain a current rotation, M15's 
initial rotation rate must probably have been even larger than observed 
today (see Einsel & Spurzem 1999). Hachisu (1978; 1982) found, using 

17 We saw a similar interesting wild-grown of theories in the gamma-ray burst community 
a few years back. 
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gaseous cluster models, that an initially rotating cluster tends to evolve 
into a 'gravo-gyro catastrophe' which drives the cluster into core collapse 
far more rapidly than would occur in a non-rotating system. If the gravo- 
gyro-driven core collapse occurred within 25 Myr, a collision runaway might 
have initiated the growth of an intermediate mass black hole in the core of 
M15. 

1.4.4 The gravitational wave signature of dense star clusters 

The peak amplitude of the gravitational wave form produced by black-hole 
inspiral is (Peters 1964) 

Here the "chirp" mass is M c hi rp = (mim2) 3 / 5 /( TO i + TO2) 1 / 5 for a binary 
with component masses toi and rri2. The frequency of the gravitational 
wave is 2/P or b, where P or b is the orbital period. The first LIGO interfer- 
ometer is expected to achieve a sensitivity h of 10 -20 to 10~ 21 at its most 
sensitive frequency around 100 Hz (Abramovici et al. 1992), corresponding 
to an orbital period of about 20 ms. The details of the wave form and the 
recovery of the signal from the noisy data complicates matters somewhat. 
With a sensitivity of h = 5 x 10 -21 , black- hole binaries will be detectable 
out to a distance of about 100 Mpc, and intermediate mass black holes to 
several kpc distances. Within a volume of §7r<i 3 we then expect a detec- 
tion rate for the first generation of interferometers of about 0.3 /iq per year. 
For ho ~ 0.72 (Freedman et al. 2001), this results in about one detection 
event per decade. LIGO-II is tentatively expected to see out to an effective 
distance about ten times farther than LIGO-I and be operational between 
2005 to 2007 (K. Thorne, private communication). This would result in a 
1000 times higher detection rate, or several events per month. 

The current best estimate of the maximum distance within which 
LIGO-I can detect an inspiral event is 

i? off = 18Mpc (^)° /6 (1-48) 

(K. Thorne, private communication). For neutron star inspiral, mi = m2 = 
1.4M , so M chirp = 1.22 M Q , Reg = 21 Mpc. For black-hole binaries with 
mi = m 2 = m hh = 10 M , we find M chirp = 8.71 M Q , R cS = 109 Mpc. 

1.4-4-1 The gravitational wave signature for intermediate mass black holes 

In S ll.4.1l we discussed a model for producing X-rays from an intermediate 
mass black hole, in this § we continue that discussion but then in relation 
to gravitational wave detectors. 
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Figure 1.23. A selected sample of binaries in the gravitational wave strain 
and frequency domain with the LISA and LIGO noise curves over-plotted. We 
assumed a distance to the source of lOkpc. The lower solid curve to the right 
is for a 2 Mq donor, the dashes for a 5 Mq donor, the dotted curves are for a 
10 M Q donor with a 100 M Q black hole (lower dotted curve) and 1000 M black 
hole (upper dotted curve). The upper solid curve is for a 15 Mq donor. Except 
for the lower dotted curve are all calculations for a 1000 Mq donor. Notice the 
enormous difference for the 2 Mq donor, which is caused by the lack of growth 
in size of the donor. 



Figur eTOSl presents the evolution of the gravitational wave signal for 
several systems which start mass transfer at zero age. Assuming a standard 
distance of lOKpc, the binaries which start mass transfer at birth emit 
gravitational waves at a strain of about log h ~ —20.2 (almost independent 
of the donor mass). The frequency of the gravitational waves is barely 
detectable for the LISA space based antennae (see Figure li~2*Sll . During 
mass transfer the binary moves out of the detectable frequency band, as 
the orbital period increases. Once the donors has turned into a white dwarf, 
as happens for the 5M© and 10 M Q donors, the emission of gravitational 
waves brings the two stars back into the relatively high frequency regime 
and it becomes detectable again for the LISA antennae. This process, 
however takes far longer than a Hubble time. 

The binary in which a 2 M main-sequence star starts to transfer 
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Figure 1.24. Evolutionary track in gravitational wave space (frequency and 
strain) for an intermediate mass black hole formed in a young and dense star 
cluster at a distance of lkpc. the initial conditions were to mimic the M82 star 
cluster MGG-11 (Portegies Zwart et al. 2004). The bullets indicate the moments 
of output one million years apart. They are connected with a dotted line. The 
intermediate mass black hole acquires a stellar mass black hole in a close binary 
near the end of the simulation at about 90 Myr which puts it in the LISA band. 



to a 1000 M Q black hole, however, remains visible as a bright source of 
gravitational waves for its entire lifetime (see Figure n~2"3*|) . The reason 
for this striking result is the curious evolution of the donor as it transfers 
mass at a slow rate. The 2 M Q main-sequence donor becomes fully mixed 
after the hydrogen fraction exceed about 66 per cent. As a result, the 
star remains rather small in size and therefore the orbital period during 
mass transfer remains roughly constant. This binary remains visible in 
the LIGO frequency regime for its entire main-sequence lifetime, only the 
gravitational wave strain drops as the donor is slowly consumed by the 
intermediate mass black hole. Mass transfer in this evolution stage is rather 
slow causing the X-ray source to be transient. Such a transient could result 
in an interesting synchronous detection of X-rays and gravitational waves. 
The star cluster contains many of binaries, some of which may have 
orbital periods small enough, and component masses large enough to be 
visible by advanced gravitational wave detectors. In figur eTOSl we show 
the population of binaries in gravitational-wave space (strain h versus fre- 
quency /) for the compact binaries for one of our MGG-11 simulations (see 
5 11.3.3(1 at an age of 100 Myr. This simulation model had initially 20% of its 
128k stars in a binary system. The majority of systems is not even on this 
graph, as they comprise of main-sequence stars, giants or other 'large' stel- 
lar object. Only binaries with at least one compact objects are presented, 
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Figure 1.25. Population of compact stars in binaries in gravitational wave space 
(frequency and strain) in a simulated star cluster at a distance of lkpc. The 
initial conditions were to mimic the M82 star cluster MGG-11 (Portegies Zwart 
et al. 2004). The cluster was 100 Myr when this image was taken. The bullets 
indicate binaries with one or two black holes, open circles indicate neutron stars 
and small dots are for white dwarfs. Where symbols overplot exactly, as is quite 
common for a circle with a small dot in the middle indicate close binaries with 
two different compact objects, in these cases they are (ns, wd) binaries. Only 
a few white dwarf binaries may be visible in the LISA band, black hole and 
neutron star binaries are generally too wide to be detectable. The intermediate 
mass black hole is not shown here, but in figure XT. 241 



and of these only the white dwarf binaries make it in the LISA band. Many 
of these systems, however, will be unobservable in realistic star clusters, as 
the star cluster is much further than the here adopted 1 kpc or the sys- 
tem ends-up in the confusion limited noise of the white dwarf population 
from the Milky-way Galaxy (see e.g. Hogan & Bender 2001; Seto 2002; 
Nelemans et al, 2004, but see Benaquista et al 2001). 
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1.4-4-2 The merger rate of stellar mass black holes 

In S ll.3.4l we discussed the evolution of relatively low density clusters through 
phase B and C. These clusters tend to eject their lack holes, in part in the 
form of binaries consisting of two black holes. The orbital parameters of 
these binaries are such that merger may occur in a relatively short time 
scale; well withing a Hubble time. We now will discuss the merger rate per 
unit volume 1Z which we predict from various types of star clusters. 

A conservative estimate of the merger rate of black- hole binaries formed 
in globular clusters is obtained by assuming that globular clusters in other 
galaxies have characteristics similar to those found in our own. Using the 
galaxy density in the local universe (see Tab. ll.2|) the result is 

TZ GC = 5.4 x l(T 8 /i 3 yr" 1 Mpc -3 . (1.49) 

Irregular galaxies, starburst galaxies, early type spirals and blue el- 
liptical galaxies all contribute to the formation of young dense clusters. 
In the absence of firm measurements of the numbers of YoDeCs in other 
galaxies, we simply use the same values of Sn as for globular clusters (see 
from Tab. ll.2(l . The space density of such clusters is then 

0YoDcC=3.5/i 3 Mpc- 3 , (1.50) 

and the black hole merger rate is 

"KyoDcC = 2.1 x l(T 8 /i 3 yr" 1 Mpc~ 3 . (1.51) 

For purposes of comparison, we also estimate the contribution to the 
total merger rate from galactic nuclei, neglecting obvious complicating fac- 
tors such as the presence of central supermassive black holes (Sage 1994), 
which may inhibit the formation, hardening, and ejection of black-hole 
binaries. The space density of galactic nuclei is only 18 about 

<j) GN = 0.012 h 3 Mpc- 3 , (1.52) 

and the corresponding contribution to the black hole merger rate is 

TZgn = 2.5 x 10~ 9 /i 3 yr" 1 Mpc~ 3 , (1.53) 

which is negligible compared to the other rates. 

Based on the assumptions outlined above, our estimated total merger 
rate per unit volume of black-hole binaries is 

K = 7.5 x l(r s h 3 yr -1 Mpc~ 3 . (1.54) 

18 Each galaxy may contain hundreds of globulars but has only one nucleus. 
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However, this may be a considerable underestimate of the true rate. First, 
as already mentioned, our assumed number (~ 10~ 4 iV) of ejected black- 
hole binaries is quite conservative. Second, the correlation between orbital 
eccentricity and binding energy and the excess of high-eccentricity binaries 
mentioned earlier both favor more rapid inspiral, causing a larger fraction 
of the black-hole binaries to merge. Third, the observed population of 
globular clusters naturally represents only those clusters that have survived 
until the present day. The study by Takahashi & Portegies Zwart (2000) 
and Baumgardt & Makino (2003) indicates that ~ 50% of globular clusters 
dissolve in the tidal field of the parent galaxy within a few billion years of 
formation. We have therefore underestimated the total number of globular 
clusters, and hence the black-hole merger rate, by about a factor of two. 
Fourth, a very substantial underestimate stems from the assumption that 
the masses and radii of present-day globular clusters are representative 
of the initial population. When estimated initial parameters (Table im 
bottom row) are used, the total merger rate increases by a further factor of 
six. Taking all these effects into account, we obtain a net black-hole merger 
rate of 

ft~6 x 10 _7 /i 3 yr" 1 Mpc -3 . (1.55) 

We note that this rate is significantly larger than the current best estimates 
of the neutron-star merger rate (e.g. Phinney 1991; Tutukov & Yungelson 
1994; Burgay et al. 2003; and many others). Since black hole mergers are 
also "visible" to much greater distances, we expect that black hole events 
will dominate the LIGO detection rate. 

1.4-4-3 Gravitational inspiral 

An approximate formula for the merger time of two stars due to the emis- 
sion of gravitational waves is given by Peters (1964): 

i -- , *'0'fe)' |i -' ,|W ' (L56) 

The sixth column of Table lTTTl lists the fraction of black-hole binaries which 
merge within a Hubble time due to gravitational radiation, assuming that 
the binary binding energies are distributed flat in log Ef, between 1000 kT 
and 10000 kT, that the eccentricities are thermal, independent of Et,, and 
that the universe is ~ 13Gyr old (Ho = 72 ± 8 kms -1 Mpc -1 ; Freedman 
et al. 2001). The final column of the table lists the contribution to the 
total black-hole merger rate from each cluster category. 

For black- hole binaries with mi = 7712 ~ 10 M Q we expect a LIGO-I 
detection rate of about 3h 3 per year. For h ~ 0.72 (Freedman et al. 2001), 
this results in about one detection event annually. LIGO-II should become 
operational by 2007, and is expected to have i? c ff about ten times greater 
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than LIGO-I, resulting in a detection rate 1000 times higher, or 2-3 events 
per day. 

Black-hole binaries ejected from galactic nuclei, the most massive glob- 
ular clusters (masses £ 5 x 10 6 M Q ), and globular clusters which experience 
core collapse soon after formation tend to be very tightly bound, and merge 
within a few million years of ejection. These mergers therefore trace the 
formation of dense stellar systems with a delay of a few Gyr (the typical 
time required to form and eject binaries), making these systems unlikely 
candidates for LIGO detections, as the majority merged long ago. This 
effect reduces the current merger rate somewhat, but more sensitive future 
gravitational wave detectors may be able to see some of these early uni- 
verse events. In fact, we estimate that the most massive globular clusters 
contribute about 50% of the total black hole merger rate. However, while 
their black-hole binaries merge promptly upon ejection, the longer relax- 
ation times of these clusters mean that binaries tend to be ejected much 
later than in lower mass systems. Consequently, we have retained these 
binaries in our final merger rate estimate (Eq. I1.55fl . 

Finally, we have assumed that the mass of a stellar black hole is 10 M Q . 
Increasing this mass to 18 M Q decreases the calculated merger rate by about 
50% — higher mass black holes tend to have wider orbits. However, the 
larger chirp mass increases the signal to noise, and the distance to which 
such a merger can be observed increases by about 60%. The detection rate 
on Earth therefore increases by about a factor of three. For 6 M Q black 
holes, the detection rate decreases by a similar factor. For black-hole bi- 
naries with component masses £ 12 M the first generation of detectors 
will be more sensitive to the merger itself than to the spiral-in phase that 
precedes it (Flanagan & Hughes 1998a, 1998b). Since the strongest sig- 
nal is expected from black- hole binaries with high-mass components, it is 
critically important to improve our understanding of the merger waveform. 
Even for lower-mass black holes (with rribh ~ 10 M©), the inspiral signal 
comes from an epoch when the black holes are so close together that the 
post-Newtonian expansions used to calculate the wave forms are unreliable. 
The wave forms of this "intermediate binary black hole regime" (Brady et 
al. 1998) are only now beginning to be explored. 



1.5 Concluding remarks 

We approach the end of the chapter on the ecology of black holes in dense 
star clusters, as part of the Como advances school in astrophysics. Here 
we will summarize a few of the key-concepts discussed in this chapter. 

Star clusters go through three (or four) evolutionary phases called A, B 
and C, each of which is dominated by either stellar mass loss or relaxation 
(see 5 11.21) . In this discussion we assumed that external influences are not 
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particularly competitive in effect, but if they are, we call it phase D. Phase 
D generally results in an early termination of the cluster. 

Black holes tend to behave differently in each of these phases. In phase 
A, it is possible to built-up a 0(1000) M© star which ultimately collapses 
to a black hole of intermediate mass (see 5 11. 3.1(1 . 

This black hole may subsequently capture a stellar companion by tidal 
effects and turn into a bright X-ray source (see 5 11.4.1(1 . 

In clusters where a runaway collision is prevented by stellar evolu- 
tion, the most massive stars collapse to stellar mass black holes in usual 
supernovae (see § 11. 2.2(1 . 

The stellar mass black holes sink to the cluster center by dynamical 
friction, pair off in binaries and eject each other from the stellar system. 
This scenario elegantly explains the absence of black hole X-ray transients 
in globular clusters (see 5 11. 3. 4|) . 

The ejected black hole binaries spiral-in due to gravitational wave ra- 
diation and ultimately coalesce. Upon coalescence they produce a strong 
burst of gravitational waves which can be detected to a distance of ~ 
100 Mpc (see §[1333 • 
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